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PREFACE 


This report summarizes the theory for a three dimensional method of moments triangular 
patch code for RCS and antenna analysis. Diagnostic capability incudes scattering image 
analysis, current display, and near field computation. The algorithm was designed to be cpu 
time and memory efficient Individual parts of this work have been documented previously 
as separate tasks. MOM3D now contains most of the original goals, thus it is time to bring 
together the appropriate theory in one report 

The FORTRAN code itself is described elsewhere. Computer codes are never static. They 
are continually changed by interested users. User manuals become dynamic documents 
compared to theory. 

This work is in support of the Advanced Aircraft Branch of the NASA Langley Research 
Center. Mr. Noel Talcott is the NASA point of contact 

This report prepared by John Shaeffer, Denmar, Inc., Atlanta Office, 3278 Hunterdon Way, 
Marietta, Georgia 30067, (404) 952-3678. 
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SECTION 1 


OVERVIEW and SUMMARY 

What Is MOM3D ? : MOM3D is a FORTRAN algorithm that solves Maxwell’s equations 
as expressed via the electric field integral equation for the electromagnetic response of open 
or closed three dimensional surfaces modeled with triangle patches. Two joined triangles 
(couples) form the vector current unknowns for the surface. Boundary conditions are for 
perfectly conducting or resistive surfaces. The impedance matrix represents the fundamental 
electromagnetic interaction of the body with itself. A variety of electromagnetic analysis 
options are possible once the impedance matrix is computed. 

What Does MOM3D Compute ? : The following analysis options are available once the 
impedance matrix is computed for each frequency and solved: 

* Backscatter radar cross section (RCS) for two orthogonal 
linear polarizations as a function of user specified azimuth and 
/ or elevation angles; 

* Bistatic radar cross section for two orthogonal linear 
polarization for user specified excitation angle and polarization. 

Bistatic RCS is computed over user specified azimuth and / or 
elevation angles; 

* Antenna pattern prediction for user specified body voltage 
excitation ports. Isotropic gain is computed for user specified 
polarization and azimuth and / or elevation angles along with 
antenna input impedance; 
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* RCS image prediction shows RCS scattering center locations 
(mechanisms) for user specified plane wave excitation of the 
body (polarization and azimuth and elevation). Either down 
range or 2D range / cross range images may be specified for 
user specified spatial extent and resolution to A./2; 

* Surface currents excited on the body as induced by specified 
plane wave excitation (polarization and azimuth and elevation); 

* Near field computation for the electric field on or near the 
body as produced by user specified incident plane wave 
(polarization and azimuth and elevation). 

MOM3D allows users the ability not only to compute RCS or antenna patterns, but also has 
diagnostic capability to aide in understanding electromagnetic responses. 

Why Is MOM3D Different ? : MOM3D is different from other surface patch codes based 
on the same basis expansion functions in the following respects:. 

* A wider choice of analysis options: Backscatter RCS, Bistatic 
RCS, and Antenna gain; 

* Diagnostic capability for: 1) Image analysis in one or two 
dimensions, 2) Surface current display, and 3) Near surface 
electric fields; 

* Efficient formulation to reduce required computer resources. 

This means that a given set of computation resources may be 
applied to electrically larger bodies. Explicitly: 1) A symmetric 
matrix formulation is used to reduce the matrix solve time by 
1/2. This varies as N 3 where N is the number of unknowns and 
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thus represents a large part of the solution time; 2) The 
symmetric matrix formulation requires 1/2 less CPU time and 
memory to compute and store the matrix elements since only 
the upper half of the matrix is used. Memory requirements vaiy 
as N 2 and become the limiting factor for electrical problem size; 

3) A user option for computing geometrically symmetric bodies. 

This reduces the matrix solution time by another factor of four 
and matrix memory storage by a factor of two ; and 4) A 
centroid approach for computing matrix interaction elements for 
all triangle to triangle distances greater than ten triangle 
dimensions. For closer spacing the user has three options for 
matrix element integration: two numerical / analytical 
approaches or the centroid approach. 

The major efficiency for applying MOM3D to electrically large bodies is for computing 
symmetric bodies. Then, compared to the patch code described in [1], MOM3D represents 
a factor of EIGHT less in matrix solution time and a factor of FOUR less in required 
memory to store matrix elements. Solution time and memory are always the limiting 
parameters for electromagnetic analysis codes. MOM3D is designed to use these resources 
in an efficient manner allowing larger bodies to be analyzed. 

This report discusses the theory behind the analysis as the reader can figure out by looking 
at the Table of Contents. 
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SECTION 2 


GENERAL METHOD OF MOMENTS THEORY 

This section discusses the fundamental theoiy for solving Maxwell’s equations as expressed 
in the electric field integral equation. This theoiy may then applied to any geometry, e.g., 
wires, 2D strips, bodies of revolution, and of course 3D triangle surface patches. References 
[2] through [7] contain more detail. The approach followed here is due to Mautz and 
Harrington [5] since they derive a very general purpose expression for matrix elements which 
then may be applied to many types of geometries and for many variations of basis expansion 
and weight functions. 

Which Integral Equation To Use? Maxwell’s four differential equations relate magnetic and 
electric fields to their electric and magnetic sources of currents and charges. These 
differential equations may be transformed to two integral equations for scattered electric and 
magnetic fields. They relate the scattered field to surface or volumetric source currents and 
charges. Green’s vector theorem is used to derive the integral form of Maxwell’s equations 
[7]. An alternative approach is to consider the scattered fields as derived from vector and 
scalar potentials. 

The usual fashion for expressing fields when dealing with scattering problems is to express 
the total field, for which boundary conditions apply, as the sum of an incident field, which 
is what exists at a given location independent of the scattering body, and a scattered field 
that is caused by the currents and charges on the scattering body. Thus 

Z To * . g 1 * + 


a™ .a*, a 


(2-1) 



The incident fields are usually plane waves as specified by propagation direction, wavelength, 
and polarization. Spherical or cylindrical wave sources could just as easily be used. E** is 
the mathematical approach of specifying the field due to outside sources. We do not know 
the source current distribution, just the resulting field. 

The integral form of Maxwell’s equations, called the Stratton Chu equations, use the surface 
parallel and perpendicular field components as source terms. These are defined as surface 
magnetic and electric currents and charges. The total surface fields are defined in terms of 
electric and magnetic surface currents and charge densities. If n is the unit vector 
perpendicular to the surface, then: 


M - E l 


A x 8; 

p* * * A • £ 

£ x A; 

p" - H x - A B 


(2-2) 


where J and M are the electric and magnetic currents while p is the corresponding charge 
density. Note that electric currents are proportional to tangential magnetic fields while 
magnetic currents are proportional to tangential electric fields. The field quantities are the 
total fields at the surface. 

The Electric Field Integral Equation (EFIE) is the Stratton Chu integral equation for the 
scattered electric field: 


= - f [jupjg - MxVg - £ Vg ] dS (2-3) 


The source terms are the electric and magnetic currents and electric charge density. 
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The free space Green’s function g relates source quantities to field quantities via phase (time 
delay) and 1/R spatial decay. It is 


8 


e -jt* 

4nR 


(2-4) 


where the wave number vector k points from source to field points and has the magnitude 
k = 2 ic / A.. The distance R is from the source point to field point location 


R - | f - 7' | * yj (x-x 1 ) 2 + (y-y) 2 + (z-z) 1 


(2-5) 


The Magnetic Field Integral Equation (MFIE) is the Stratton Chu integral equation for the 
scattered magnetic field: 


ff*V,) - f I -y'ueMg + JxVg + -H^Vg ] dS (2-6) 


The source terms are the magnetic and electric currents and magnetic charge density p m . 

For perfect electrically conducting bodies the magnetic currents are identically zero since 
there can be no tangential electric field. Then either the EFIE or MFIE may be used. 
Typically the EFIE approach is used for open and / or closed surfaces while the MFIE is 
used only for closed surfaces (due to problems in numerically evaluating the Green’s 
function). 

MOM3D uses the EFIE approach with electric current unknowns. 
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Magnetic currents must be included for the more general case when considering magnetic 
materials, e.g., radar absorbing material. This doubles the number of unknowns, with matrix 
element memory increasing by four N 2 , and matrix solution time increasing by eight N 3 . We 
also must formulate the system matrix to include the magnetic sources. This requires 
additional analytical and computation effort 

For Perfect Electric Conductors, PEC, there are no magnetic source currents or charges 
since the tangential E field is zero. This greatly simplifies the EFIE and MFIE because now 
we only have to solve for the unknown surface currents J. 

The astute reader will ask why we also do not explicitly solve for the charge density. The 
reason is that continuity relates current to surface charge density, i.e., when current flows 
into a region the charge density increases while when current flows out of a region the 
charge density decreases. Mathematically, 


V-7- = -/«p (2-7) 

dt 


where we have explicitly assumed time variation of e i “‘, i.e., we assume a time harmonic 
monochromatic frequency. 

Boundary Conditions: The next major step is to specify the boundary conditions on the 
tangential electric field. For resistive boundary conditions the tangential electric field is 
equal to the surface impedance Z** times the current density J. If the surface is a perfect 
conductor, - 0, then we have the requirement that the tangential E field vanish. When 
the field observation point is on the body surface we have 

* X - A X ( £* * £•“ ) - z^J (2-*) 
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The known quantities are the incident field and surface impedance. The unknown is the 
current density J (which also appears under the integral sign of the EFIE). The boundary 
conditions lead to a Friedholm integral equation. When Z*, rf is zero, this is an integral 
equation of the first kind, ie., the unknown J appears only under the integral, or when non 
zero, an integral equation of the second kind since J also appears outside the integral. 

It is customary to rearrange the unknowns to the left and known quantities to the right: 

-j! x -Axi* ( 2 * 9 > 


It is also customary to write the integral expression for E* 3 * as a linear operator L operating 
on the current J, as L(J): 


V •/ 


IXj) - -£~ =;*T| / fjg - -fr ] dS 


( 2 - 10 ) 


where kq = cop, toe = k/r| and have set the magnetic sources to zero. 

Matrix Solution of the Integral Equation: The formal solution of the integral equation 
yields a matrix or set of simultaneous equations. Sub domain basis functions with unknown 
coefficients j represent the unknown surface current: 

/ ■ E ho < 2 - n > 

l-I 
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The sub-domain basis function f(t) contains the vector and functional characteristics for the 
local surface current expansion function. This representation of a function is a 
generalization of other types of series expansions such as Taylor or Fourier. 

The boundary conditions are enforced at each point or region on the body to determine the 
unknown currents. The boundary conditions, where we know the solution, are enforced at 
N unique locations on the body, i.e., the tangential field must be equal to Z** J. At each 
location on the body the total field is the sum of the incident field plus the field caused 
(scattered) by currents located on every other part of the body. This leads to a coupled 
system of equations by which every part of the body affects every other part ( Maxwell’s 
equations are elliptic ). 

Impedance matrix elements express mutual body coupling. When enforcing the BC at 
specified locations on the body, we sample the region surrounding the BC location. 
Mathematically we obtain an average value by using a weight function for each BC 
enforcement location. Explicitly we multiply the boundary condition equation by a surface 
tangential vector weight function W and use an inner product integral defined as 


- / *, • dS p 


( 2 - 12 ) 


where the integration is over the p th surface area which surrounds the p th BC match point. 

The current series expansion is inserted into EFIE expression for scattered field and the 
boundary conditions applied to obtain the matrix form of the integral equation. 

The inner product is a linear operation which slides over the integrals. A set of N equations 
for N unknown current coefficients results: 


2-6 



£ + Z^<W p J>j, - <#,,£*"> for p - 1 to N (2*13) 

«-i 

for N unique locations on the body. A matrix, Le n set of simultaneous equations, represents 
the original elliptic integral equation. 

The matrix element interaction between body surface patch areas p and q is the inner 
product < Wp, L(fq) > = Zpq. The voltage vector solution forcing function is the right hand 
side inner product, V q = < Wq, E** >. 

With this notation, the original integral equation becomes a matrix equation: 


z n z i2 ••• z in 


it 


v i 

• • • 


h 

= 

V 2 

• • • 
ha ••• z nn 


J'n 


V H 


(2-14) 


or in compact matrix notation 


ZJ = V 


(2-15) 


where Z is a N by N square matrix, J is the unknown current vector, and V is the voltage 
vector forcing function. J and V are column vectors with N rows. 

The solution is obtained by computing the matrix elements Zp, and solving for the currents 
J for various body excitations forcing functions V, e.g., plane waves or antenna voltage ports: 
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(2*16) 


Voltage Vector: The forcing function for scattering problems is usually specified as an 
incident plane wave. The tangential component of the incident field is automatically 
obtained from the surface vector weight function W which is defined tangent to the body 
surface. The excitation for RCS computations is usually a plane wave: 

E* ■ E^u'e-*-* <2-17) 


where Eg is the wave amplitude, u is the polarization unit vector, k"' is the vector direction 
of the plane wave, and R is any coordinate location vector. Planes of constant phase are 
perpendicular to k“. The components of V for plane wave incidence are then: 


v‘ = < w q , > = f ( w t • a m ) dS q 


(2-18) 


where u" “ 8 or * is the polarization unit vector, k" is the wave number in the direction of 
propagation of the incident wave, and R, the local coordinates of the q th sub-domain. 

Matrix Elements: Mautz and Harrington, [5], derived a universal formula for computing 
matrix elements. This recipe is applicable to any geometiy for arbitrary choice of expansion 
functions f and weight functions W. The expression for matrix elements is: 
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Z„ - <#, uj) » ■ <7*11 //[<#, •/,) * #, v* 1 as, as, <*•»> 


which for three dimensional surfaces becomes a four fold integral, i.e., surface to surface. 
Reference [5] used the 2D surface divergence theorem on a closed surface to convert the 
second term dot product: 


f Vg dS - -f g V#dS 


( 2 - 20 ) 


to arrive at the universal matrix element interaction expression: 


Z = +jkr\ f dS / dS [ W •/ " — < 2 ’ 21 > 


This expression is the "mother of all matrix element formulas." It applies to all geometries 
and to all types of expansion functions f and weight functions W. It applies to bodies of 
revolution (BOR’s), wires, 2-d strips, 3-D patch, wires attached to BOR’s or patches, BOR’s 
coupled with 3-D patches, etc.. 

Scattered Field: The heart of the MOM solution is to solve for the current distribution on 
the body. From these currents one is then often interested in computing the resulting 
radiation pattern, E**, from which we obtain radar cross section or antenna gain patterns. 
Far field radiation is given by the transverse components of the electromagnetic vector 
potential A: 
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( 2 - 22 ) 


£«* JM '(0,(|>) - -ja(A'X) = -y«|i f (n«J) g dS 



Figure 2-1 Far field distance approximation geometry 

where a = 0 or <p polarization. The standard approximation for the distance from source 
to field point is then made using the law of Cosines, Figure 2-1, which is valid whenever the 
field point distance is much greater than the local body dimensions: 


R * I* -*,l - *0 


, - * 


2 ' 

*0 




11/2 


(2-23) 


* Ro ~ Ro 'Rp 


where Ro is the distance to the far field point, Rp is a local source point on the body, and 
Ro points in the same direction as k sat . The distance R in the Green’s function is 
approximated by Ro in the denominator (amplitude) and by the two term approximation for 
the phase term. The series current basis function expansion is used for the current J. The 
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final expression for scattered field then becomes a complex dot product of the row 
measurement vector R and the solution column vector j: 

t 


£‘•**(6, 40 = 3^- e' Jk *• R a • J 


(2-24) 


where the row vector R is 


J£( 6,<t>) « / ( A* ‘7 e* s£ ~*> dS p (2-2 5 ) 


where the f p is the current expansion function for the p th sub-domain, n* " 6 or * is the 
polarization, k 5 ”* is in the direction of the scattered field, and Rp is the local coordinates of 
the p th sub-domain. Harrington [6] derives the same formulation using reciprocity. 

For backscatter k** = - k” therefore the expression for scattered field is the same as the 
voltage vector expression. If f = W, the same computation routine may be used to compute 
either. 

Computer Resource Considerations: The number of unknowns in a moment solution is 
proportional to the body size as measured in wavelengths since the currents must be sampled 
on a wavelength scale, typically 7 to 10 per wavelength. Computer resources of memory and 
available cpu time limits the body size which one may compute. 

Surface currents require vector unknowns. Each sub-domain area requires two scalar 
unknowns, one for each basis vector direction. Let us characterize a body by dimension L 
and area L 2 . If we require n unknowns per linear dimension, then a surface with area A 
would require 2n 2 unknowns per square-wave length (98 to 200 unknowns) of area. Matrix 
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storage then would be proportional to ( (n L) 2 ) 2 = n 4 L 4 . Matrix solution time, which 
varies as N 3, would be proportional to ( (nL) 2 ) 3 = n 6 L 6 . Thus as the body grows in 
electrical size, required memory grows as (linear size) 4 and matrix solution time grows as 
(linear size) 6 . 

Linear sample density n should be as small as practical and still yield a satisfactory solution. 
This will reduce memory requirements (proportional to n 4 ) and matrix solution time 
(proportional to n 6 ). Computer resource requirements grow very rapidly with linear body 
size L and sample density n. We thus see the importance of optimizing the algorithm to 
conserve computer memory and time. 

It should be noted that a fundamental reformulation of the analysis will be required if we 
are to compute very large electrical bodies. Alternately parallel computer architecture will 
be required. 

The impedance matrix is of rank N 2 , i.e., N unknowns. This means that at least N 2 memory 
locations be available to store the complex matrix elements. This is the dominant 
requirement for memory storage. Matrix fill cpu time will also be proportional to N 2 . 

Matrix solution time is proportional to N 3 . Matrix inversion (or equivalent process) 
dominates the solution time for large problems. Various approaches can change the 
proportionality constant (as optimized in MOM3D). For a given geometry and frequency, 
the impedance matrix may be computed, solved, and then saved to disk storage. Then 
various analysis options may later be performed using the saved LU decomposed matrix. 

Matrix solution may be obtained in principle by a variety of approaches: inverse, Gaussian 
elimination, Crammer’s rule, LU decomposition, iterative solutions, etc.. As discussed in [2], 
[8], and [9], LU matrix factorization requires 1/3 the time of direct matrix inversion. 
Therefore MOM3D utilized the LU factorization approach. 
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The LINPAK [9] library of matrix routines for LU decomposition and back solution was 
used. The matrix formulation was specifically made symmetric by choosing Galerkin weight 
functions, i.e., f = W and the LINPAK library routines for symmetric matrices was used. 
This required 1/2 less matrix storage and computation time, and, more important, reduced 
the LU decomposition time by 1/2, which for large problems is significant. 

Future adaptations could take advantage of parallel processor architectures using matrix 
libraries specific to these machines, such as the LAPACK library now under development, 
reference [10]. 

When using LU decomposition, the back solution time to compute the currents is 
proportional to N 2 . For backscatter RCS problems, a new voltage vector and current 
solution is required for each excitation angle. Therefore the solution time becomes 
proportional to M_ c ,_ N 2 . Since we typically compute backscatter every degree, M may 
equal 181 to 361 depending on the angle range. Thus for small bodies, if N is less than M, 
the back solution time may become greater than the LU decomposition time. 

The time to compute the voltage or row vectors is proportional to the number of unknowns 
N. This is not a significant time requirement compared to back solution time. Other 
computation overhead is required to compute the required geometry parameters from the 
input geometry description. 
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SECTION 3 


MOM3D SURFACE PATCH FORMULATION 

The theory developed in Section 2 is completely general. This section will apply the theory 
to the surface triangle basis functions introduced by Rao, Wilton, and Glisson [1]. Matrix 
and body symmetries are used in MOM3D to create a computational efficient code 
compared to [1]. MOM3D has diagnostic capabilities for imaging, near fields, and 
currents, as well as for backscatter & bistatic RCS and antenna gain pattern computations. 

Surface patch basis functions must: 1) model two independent surface vector directions; 2) 
have desirable derivative properties since the second term in the matrix element formula 
requires the divergence of J; and 3) model arbitrary curved or planar surfaces. A 


h*lMt.g«o n - 115.0 EL - 30.0 0 [STANCE - 100.0 MIRRORED 



HELMET j d • 2.04'. • - 6.12* Bc„ • 1.0 

REFLECTION AXIS : 1 

NuatMT of Point*. NPOINT.i 145 
■ o-f Trlftngl««, NLTRIA-t 262 
Unit*, l»M. 3* Inch**: 3 

Sur-ftc* lnp«dtne« (R*. In) t 000.0 00.00 

Figure 3*1 Arbitrary surface modeled by triangular patches 
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rectangular surface patch basis function code is reported in [11] but it has obvious problems 
when modeling arbitrary curved surfaces. These three criteria were used by [1] to develop 
the surface triangle basis functions. They can model a wide range of geometries as 
illustrated in Figure 3-1. Both open and closed geometries can be modeled. 

Current Couple: A pair of coupled triangles 
that share a common edge were introduced 
by [1] as the fundamental vector basis 
function. Figure 3-2. The vector current 
flow is from one vertex flowing across the 
common edge toward the opposite vertex. 

Each triangle pair is a couple and is the 
fundamental unknown quantity in the 
analysis. This basis function is such that its Figure 3-2 Triangle couple basis function 
magnitude is like a triangular hat function. 

It has the value zero when r is at the vertex, unity when r is anywhere on the common edge, 
and zero when r is at the opposite vertex. The vector basis function associated with the q 
th edge is then 
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By arbitrary definition, current flow is from the first triangle of area A* to the second 
triangle of area A*. The length of the common edge is 1,, p% is a vector from either vertex 
toward the common edge, and S is a unity sign showing current flow direction. This 
definition is slightly different from that in [1] in that S = ± 1 is explicitly used. This 
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simplifies the following analysis and resulting computer code since the vector p can always 
have a common definition regardless of whether a triangle is + or This basis function is 
the surface version of triangle basis functions used for wire or 2D analysis. It has the value 
of zero when p is at either vertex and a value of unity when p is on the common edge. The 
slope or derivative (divergence) is positive constant in the + triangle and negative constant 
in the - triangle. This derivative feature is such that the triangular hat current yields a pulse 
doublet charge distribution. 

The surface divergence operator for triangular patch geometry is 

A 9(p 4 ) (3-2) 

P dp 



which results in the following expression for the surface divergence of the vector basis 
function (charge density) 



(3-3) 


These vector basis functions are used to represent (approximate) the surface currents. 
Since a triangle has three sides, and if each side has a common neighbor, then it contributes 
up to three basis functions (couples). The current in the q th patch is represented by up to 
three basis vectors. These are not orthogonal as is the case of usual basis vectors, however, 
they are independent and "span the space" to represent any physical direction of current 
required by the solution of Maxwell’s equations. 
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Couple Basis Function Average (Moment): The fundamental current vector unknown is 
between two adjacent common triangles and not a single triangle by itself. It is convenient 
to describe the location of the couple centroid. The average or moment over the two 
triangles of the couple basis function is [1]: 
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where p 6 ’* is the vector from the free vertex 
to the centroid and i** is the position 
vector of the triangle centroid, Figure 3-3. 
We see that the average of the basis 
function over the triangle couple pair is 
equal to the common edge length times the 
vector between the centroid of each 
triangle. The center of this vector is the 
center of the vector basis function I”'”®*. 



Figure 3-3 Net couple moment (current) is 
from triangle centroid to centroid 


Matrix Elements: The development of 

expressions for matrix element interaction requires one to choose a specific weight function 
W and to consider if a symmetric matrix is to be obtained, i.e., Zp, = Z v . Various options 
are possible as discussed by [2] and [7]. The general expression of matrix elements for three 
dimensional surfaces, Section 2, shows that is a four fold integral, i.e., the E field 
integrated over surface patch p (with averaging weight function W) as caused by currents 
located on surface patch q with basis function f. Reference [1] simplified the computation 
by choosing delta or point matching weight functions so that a two fold integral resulted 


3-4 


(point to surface), however a symmetric matrix did not result To obtain matrix symmetry, 
the weight functions must be of the same form as the current expansion functions, W q = f q , 
which is called the Galerkin choice for W. The penalty paid to achieve matrix symmetry 
(with resultant reduction in LU decomposition and memory storage) is the requirement for 
four fold integration. As we shall shortly see, the centroid integration approximation 
eliminates this integration requirement (We can have our cake and eat it too?) 


The reader can easily see from Section 2 
matrix element expression that a Galerkin 
choice for weights yields a symmetric 
matrix. This will be our choice. All we 
have to do is to insert the appropriate 
values to obtain the required expression. In 
doing this we must remember that the 
integrations are from one pair (couple) of 
triangles to another pair (couple). Figure 3- 
4. The result is: 



Figure 3-4 Matrix element interaction is 
between pairs of triangles (four terms) 
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Since each couple has two triangles, each matrix element has four terms, i.e., each triangle 
interacting with the others, (+,+), (+,-), (-,+), and (-,-) : 
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The next step is to pull out of the indicated integration the constant parameters and to 
define the major Green’s function integrals. We define G (1) w and G (#) „ normalized to 
triangle patch areas: 
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We then arrive at the required expression for the matrix elements by using the appropriate 
expressions for f and divergence f: 


z„ - m i t i. 


s r .s t . [ 1 H G™. 

S, S r [ 1/4 G^V _ 
S,-S t . 1 1/4 G™. 
[ 1/4 G™. 


- l/* ! G® f . ] . 

- 1/t 2 G® . 1 * 

- 1/Jfc 1 G™ f . ] - 

- I /* 2 G® f . 1 ) 


(3-8) 


This is our exact (numerically) desired expression. No approximations are made at this point 
other than the choice for the current expansion functions f and weights W = f. 

Green’s Function Computation and Approximations: The principal task is now to compute 
G (1) and G (0) that integrate the free space Greens’s functions over each pair of triangles. 
This of course is a four fold surface to surface interaction. G (1) and G (0) are related to the 
electromagnetic vector and scalar potentials. We start by estimating the dominate term in 
the matrix element computation. We assume that the surface is sampled such that p < A/10 
so that the ratio of the first term in Z (vector potential) to the second term (scalar potential) 
is approximately: 


Vector = °- 2 ? L£ < *L . o.i (3-») 

Scalar l/* 2 100 

which shows that the dominant interaction between surface patches is from the scalar 
electrostatic term. Based on this observation, we approximate the G (1) integral as: 
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where p e por<I is the vector from the free vertex to the triangle centroid, i.e., we have used a 
centroid approximation for the vector dot product portion of the integrand in G (1) . 

Centroid Approximation: The computation of G w remains. This is still a surface to surface 
four fold integration involving the free space Green’s function as the integrand. The 
appendix discusses a combination analytical / numerical integration approach, i.e., one 
surface integration is done analytically and the second surface integration is done with either 
a three or one point numerical integration. This section discusses the centroid approach. 


The integrand has a 1/R intensity 
decay and phase (time) delay. 

The surface patches have linear 
dimensions less than 1/10. 

Therefore the stationary part of 
the phase involves the distance 
between triangle centroids, 

Figure 3-5, while the varying 
portion has at most a 2 (1/20) = 

A/10 = 36° variation in phase. If 

we recall the typical Figure 3-5 Distance between triangles using centroid 
approximation used in plane 

wave EM analysis, a 22.5° variation is considered constant. Thus the phase part of the 
integrand is relatively constant centered on the centroid value. 



The 1/R amplitude portion of the integrand is dominated by the distance between triangle 
centroids that is just the leading term of a Taylor scries (expanded about the centroid value). 
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Consider Figure 3-5 where the distance between source and field points, R„ is expressed in 
terms of the vector from triangle centroids, R„ e and two local vectors in each triangle, r p and 
r q . Then we have 
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With this centroid approximation, G (0) becomes a veiy simple expression, i.e., just the free 
space Green’s function evaluated at the centroid distance: 
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Notice that no integrations are required. 


Self Term Approximation: Self term matrix elements represent the field at a couple due to 
its own currents, Zpp. Self term quantities also have four terms, two self triangles, (+,+) and 
(-,-) and two adjacent triangles, (+,-) and (-,+). When evaluating the self triangle, the 
centroid distance is zero and clearly an alternative expression is required to prevent the non- 
physical case of division by zero. The approach we will use is that by Harrington [6] who 
simply integrated the free space Green’s function over a patch that he took to be circular 
with a radius defined such that the circular patch had the equivalent area to the actual 
(triangle) patch under consideration. The main contribution to the integral is from the 
region surrounding the centroid thus the result does not change appreciably provided the 
triangle aspect ratio is not extreme, i.e., the height to width ratio remains near unity. The 
self term approximation starts by writing the self term as: 
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(3-13) 
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The integrand is then expanded about the center of the triangle, R = 0, and only the first 
two terms retained (remember that R << A.): 
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The self term then becomes 
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This self term approximation along with the centroid approximation is used whenever 
MOM3D users specify centroid integration for all couples. The analytical / numerical 
approaches discussed in the appendix yield approximately the same result for self terms. 

Voltage Vector The right hand side forcing function of the matrix equation is the voltage 
vector. When computing either backscatter or bistatic radar cross section the voltage vector 
is due to the specified incident plane wave: 
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E* - &•••<** E 0 e-t** M 

where the polarization reference is to a spherical coordinate system and may be either 6 or 
the direction of propagation is k**, and Eq is the amplitude taken as either 1 volt / meter 
or in MOM3D as tj = 377 volt / meter to reference currents to a unit incident magnetic 
field. 

The voltage vector for our Galerkin weight functions, W * f, is then 
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The centroid approximation is used to compute the two fold surface integral: 


v , = k ( s; (iv-ii') ^ * s; (?;. «■) \ e 0 


(3-18) 


In MOM3D, each row of the column voltage vector is computed by this expression. 

Another approach for evaluating the voltage vector is to use f "* for the integration over the 
two triangles forming the couple: 
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where the specific form for f*** is used and the incident field phase is evaluated at the 
midpoint between triangle centroids. 

Scattered Field: The field radiated by the currents on the body is computed using the 
general expression developed in Section 2. The radiation integral is evaluated as a complex 
dot product of a row measurement vector with the column current vector. The elements of 
the row measurement vector, when Galerkin weight functions are used, W = f, has 
identically the same form as the expression for the voltage vector. The astute reader will 
notice that the elements of the row vector are defined for If* pointing away from the local 
origin while the voltage vector is elements are defined for If* pointing toward the local 
origin. This change of sign is what makes the row and voltage vectors have equivalent form. 

The reader is cautioned that while the voltage vector subroutine may be used to compute 
the row vector, V p is not equal to R p unless we are computing the backscatter case when If* 
= - If*. The elements of k, as expressed in rectangular coordinates in terms of the polar 
angle 0 and azimuth angle q> for an incident wave toward the origin, are 
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This is a radial vector pointed toward the origin. The two polarization unit vectors are 
orthogonal to k. The usual choice is to pick the 0 and <p directions, tangent to a sphere, as 


ff* = ( cos 6 cos<|> , cos0 sin4> , -sin© ) 
■ ( -sin4> , cos<(> , 0 ) 


(3-21) 
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Figure 3-6 Polarization unit vectors and 
incident field unit vector 
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SECTION 4 


GEOMETRY SYMMETRY CONSIDERATIONS 

Considerable reduction in computer resources can occur if the geometry has symmetry (see 
Miller [2]). Examples of body symmetry would be bodies of revolution, bodies with planes 
of symmetry such as aircraft, and bodies such as hexagons with multiple symmetry planes. 
Body symmetry reduces solution resources because: 

* The impedance matrix represents the interaction of each part of the body with 
another, therefore symmetry implies that the coupling between similar body 
parts is equivalent. This means we do not need to recompute those elements 
of the impedance matrix. This implies a reduction in matrix element storage 
and compute time, an N 2 dependence. 

* The solution of the impedance matrix can take advantage of equivalent sub 
matrices to reduce the solution time N 3 dependence for LU decomposition. 

Body symmetry does not automatically reduce computer resources. The solution algorithm must 
be specifically designed to incorporate symmetry. 

Body symmetry is a completely different concept than matrix symmetry obtained by rising 
Galerkin weight functions. 

MOM3D solves body geometries with either no symmetry or with a plane of mirror symmetry 
about the X, Y, or Z axis. 

With body symmetry the impedance matrix splits into two smaller matrices, each with 
approximate rank N/2 where N is the total number of unknowns on the body. The following 
efficiencies result: 


d 


4-1 


* Matrix storage reduced by two. We have two matrices of approximate rank 
N/2, therefore matrix fill time and memory storage is: 2 (N/2) 2 = N 2 / 2. 

* Matrix LU decomposition solution time reduced by four. Two matrices of 
rank N/2 require a solution time of: 2 (N/2) 3 = N 3 / 4. 

The savings in computer resources becomes very significant for large electrical bodies when 
N is large. 

Mirror Geometry Symmetry: MOM3D has 
incorporated mirror symmetry about either 
the x, y, or z axes if specified by the user. 

The geometry is described somewhat 
arbitrary as left side, right side, spine, and 
reflection plane, Figure 4-1. The left and 
right sides are those couples that are all 
either on the left or right. The spine 
couple unknowns bridge between left and Figure 4-1 Geometry example that has 
right sides, i.e, they have one triangle on s 5' mmctI ’' • bou ‘ the x " 0 (Y ‘?> P ,ane - 
the left and the other on the right side. 

The reflection plane couple unknowns are not reflected since they are on the symmetry 
plane. 

We make use of specific relationships between these four type of unknowns to express the 
system impedance matrix as two separate matrices called Z even and Z odd. 

System Matrix: The system impedance matrix is written in terms of these four types of 
unknowns: Z^ and Zrr the left and right side unknowns; Z^ the spine unknowns; and Zw 
the reflection plane unknowns. With this breakdown the system matrix equation becomes 
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Matrix element relationships: Geometry 
symmetry means that matrix elements are 

UfcCoupfe 

the same for left-left and right-right sides, rr— 
i.e., the interaction of the left side of the ^ 

body with itself is the same as the right side 
of the body with itself, = Zrr. The 
reflected couples are defined to have an 
opposite direction, Figure 4-2, so that the 
interaction of the spine with the left side is Figure 4-2 
the same as the interaction of the spine °PP^ te direction so ^ at ^ls = ^rs and A.v 

with the right side, = Z^ and the 

interaction of the left side with the reflection plane couples is the negative of the interaction 
of the right side, Zl V = - Zr V . The interaction of the spine couples, which are perpendicular 
to the reflection plane couples, is zero, 2^ v — Zys = 0. Due to matrix symmetry the 
following are equivalent: = Z^ Z^ = Z^ , Z^ = Z^ = Z^, and Zr V = Zvr. 

The proof of these assertions is left to the reader. 
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We develop the even and odd impedance matrices by multiplying out the system matrix 
equation, adding and subtracting various rows, and then use the sub-matrix equivalences. 
Multiplying out yields: 
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Wc define even and odd as follows: matrices and unknown current vectors 
and arid voltage vectors V m and V^. 

We then start by adding and then subtracting (1) and (2) to obtain: 

( Z u + Z„ )( I L * /, ) + 2 Z U I S = ( V L + V R ) (5) 

(4-3) 

( Zu - z a )( l t - /, ) + 2Z 1K / K - ( r 4 - K, ) (6) 


The even system of equations is defined as the set (5) and the spine set (3) with a factor of 
2 incorporated to simplify later equations: 
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The odd system of equations is defined as the set (6) and the reflection plane set (4) along 
with the factor of 2: 
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If the original system matrix is symmetric due to Galerkin choice of weight functions, then the 
even and odd matrices are also symmetric. 

The explicit even matrix and column vectors are: 
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while the odd matrix and column vectors are: 
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M0M3D computes Z and V in normal fashion to create Z even/odd which are then solved 
via LU matrix decomposition. The even/odd voltage vector is computed from the normal 
V and the resulting even/odd currents are computed. Finally, the left and right side currents 
are obtained for scattered field analysis from the even/odd currents: 


« c Ifmt ” I odd ) / 2 
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For bodies that have a plane of symmetry, this algorithm approximately reduces matrix 
solution time by four and matrix fill and storage requirements by two. 

Body symmetry and matrix symmetry form the cornerstone concepts for making MOM3D a 
computation efficient algorithm. 
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SECTION 5 


BISTATIC IMAGE ANALYSIS USING k SPACE CONCEPTS 

The experimental development of microwave images as pioneered by Dean Mensa [12] is 
one of the most powerful tools used to understand the scattering from various geometries. 
Imaging may be in one dimension, i.e., down range; or in two dimensions, i.e., down and 
cross range. This capability allows one to understand the scattering process in terms of 
specific scattering centers and mechanisms. Image development has been mostly 
experimental. While one could apply the same methods to predictive scattering algorithms, 
the computation burden has always been considered to great This occurs because, 
experimentally, down range information is obtained by illuminating the target over a 
bandwidth of frequencies typically numbering 16, 32, 64, 128 or even 512. To do this with 
a method of moments analysis, one would have to recompute and solve the system matrix 
for each frequency. This computation burden is so great for large problems that the swept 
frequency approach is seldom pursued. 

B. A. Cooper developed a new approach that requires only one computation of induced 
currents and therefore only one MOM matrix computation for down range images. A 
formal bistatic k space image theory was then developed by the author. This formulation 
computed cross range images without smear. The bistatic k space analytical image technique 
does not require a MOM code matrix solution for each frequency. Only one current 
distribution (matrix computation) is computed at the frequency of interest The image is the 
Fourier transform of the k space bistatic scattered radiation for values of k s “* that 
correspond to downrange and cross range. The natural Fourier transform variables are wave 
number k and spatial position R. If we compute a bistatic field as a function of k** = 
(k<k*mw*e, mge) ^ Fourier transform of the scattered field is naturally a function 
of the transformed spatial coordinates, R"* = (R dw ’ w * r , R** 0 '' ""•*). The computation of the 
scattered field in k space is a generalization of the standard bistatic radiation integral. The 
difference is that E** is computed in term of k** for down and/or down/cross range rather 
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than in terms of the usual bistatic angles (6,<|>). Body currents are computed only once at 
the user specified incident angle and polarization. 

Experimental Image Concepts: We start the discussion by looking at how experimentalist’s 
obtain microwave images. This will form a background for the k space analytical image 
development. 

The fundamental requirement for imaging is to obtain a scattering response that is a 
function of body location. This is done by causing the relative phase to change in both down 
and cross range. Down range phase is accomplished by sweeping the frequency which 
changes the relative downrange position (phase) of the scattering centers. Cross range phase 
is accomplished by rotating the body. In an electromagnetic sense, we rubber band or 
stretch the body in phase (time delay) so that we can reconstruct the physical scattering 
locations via the Fourier transform. Experimentally, the only way to move scattering centers 
in down range is to vary the number of wavelengths in a down range direction, i.e., change 
the frequency. Greater detail may be found in [12] and [13] and the author’s section of 
Chapter 8 [18] taken from RADAR REFLECTIVITY. 2nd edition. 

The limitations of experimental imaging are resolution and image smear or focus. 
Resolution increases with bandwidth: 


Ar _ c X 1 
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Typical experimental systems are usually bandwidth limited although a few experimentalists 
have pieced together body responses over several waveguide bands. A 10% bandwidth yields 
approximately a 5 wavelength resolution. At X band this is usually satisfactory, while at L 
band or VHF this clearly is not very practical. 
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Cross range resolution is increased by increasing the angular extent over which the model 
is rotated. The small angle approximation, cos 0 * 1 and sin 0*6, however, limits 
resolution. Large angle extent causes the image to smear and become de-focused. It is 
possible to refocus the image with an appropriate algorithm, but this becomes a difficult 
image processing factor. 

Additional concerns with experimental images are: 1) at what frequency does the image 
represent (since the image must be obtained over a bandwidth of frequencies)?; 2) what 
angle does the image represent for cross range since the model must be rotated through a 
finite angular extent?; and 3) what physical scattering mechanisms are we not observing or 
interpreting incorrectly due to the frequency and angular rotation sweep? 

Bistatic k Space Image Background: The initial concept was explored by Cooper who 
developed the technique for down range images for a body of revolution MOM code using 
the experimental image approach with a synthetic frequency sweep in the radiation integral. 
It was soon apparent that high resolution images in down range were possible leading to 
spatial resolutions approaching XU and most important, requiring only one MOM code 
matrix computation and inversion to compute the body currents. Very useful diagnostic 
information could be analytically computed with any MOM code for a very small additional 
computation cost. 

When this technique was applied by the author to compute two dimensional images using 
the synthetic approach with frequency sweep and body rotation as described in [18], the wide 
bandwidth and angular rotation extent violated the small angle approximation. The result 
was a de-focusing (smear) of the 2D images and scattering amplitudes. Clearly one could 
have restricted the technique to small angular sweeps; but this would have negated the 
desirable high resolution. 

To achieve high resolution in both range and cross range without image de-focus, the author 
reformulated the image mathematics from the experimental approach where frequency and 
angle are the primary variables to an analytical bistatic k space approach where the primary 
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variables are the bistatic wave number in down and cross range directions and the 
corresponding transform spatial position variables in down and cross range. 



Figure 4-5 Bistatic image k space 

A natural bistatic k vector for down and cross range corresponds to the spherical coordinates 
of ( k„ k + ), Figure 5-1. The downrange direction is radial, lq, while the cross range 
direction is in a plane perpendicular to the radial down range direction, i.e., either k, or k v 
A full three dimensional image could be computed using this approach. The present 
implementation is limited to either downrange or down / cross range. 

The bistatic k space image technique has the potential for the following (not all are 
presently implemented in MOM3D): 

* Resolution up to Ul unlimited by the usual experimental 
frequency and angle extent bandwidth concerns; 
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* Image focus / smear does not occur due to the formulation of 
the approach; 

* Images are computed at the frequency and angle of body 
excitation. The body currents are computed only for this k iDe . 

In contrast to the experimental approach, the currents do not 
change with changes in the bistatic k** vector sweep; 

* One, two, or three dimensional images may be obtained. The 
limiting feature for obtaining 3D images is the display of the 
solution since multidimensional FFT algorithms are available; 

* The scattering body is imaged in a bistatic sense from the same 
direction as the excitation, i.e., a backscatter image. However, 
a more general bistatic approach is entirely possible since the 
center of k** is not required to be the negative of k“; 

* A co-polarized image is implemented, but a cross polarized 
image could easily be computed. 

The above list of possibilities shows that the bistatic analytical approach to imaging can 
potentially yield substantially more information than experimental images. The bistatic image 
approach can be applied to any predictive algorithm for electromagnetic scattering or antenna 
radiation, e.g., Physical Optics. 

We can expect to see similarities and differences for a given target as imaged experimentally 
with swept frequency versus the analytical bistatic image. The analytical image is obtained 
from body currents that are excited only at one frequency and one excitation angle while the 
experimental image is a response over a bandwidth of frequencies and angle sector 
illumination. Experience must show how the two approaches will differ. One could argue 
that the bistatic image is a truer representation of the physical scattering mechanisms since 


5-5 



the body is truly excited at one frequency and angle. Thus line source radiation mechanisms 
such as reflected edge or surface traveling waves should image as radiative sources over the 
entire edge or surface as compared to the experimental approach where these mechanisms 
always show as emanating from an aft edge or tip vertex ( because the swept frequency 
stationary phase point here is only at the fixed geometry end region ). 

Bistatic Image Theory; The key breakthrough in developing this approach for multiple 
dimensions was to abandon the experimental frequency and body rotation approach. The 
new concept falls naturally out of standard FFT theory [14]. We start by writing the bistatic 
scattered field as a function of three spatial coordinates that naturally constitute an image, 
i.e., down range and in a plane perpendicular to down range: 


£ “ ,P ( r r» r e» r *) 
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This is the Fourier transform for the bistatic scattered field computed in k““ space, i.e., the 
square magnitude of this quantity converted to RCS is the image, i.e., intensity as a function 
of spatial position. 

The variables are: 

| E*’ p (r r ,r # ,r^ )| 2 is proportional to the amplitude of the desired image, e.g., 

RCS; 

k”* = k (do,raraD * c ’ cro " “s** ss (kpk^k^) is the wave number in spherical 
coordinates corresponding to down range and / or the two cross range 
directions; 

E^k^k**) is the scattered field in If* space computed from the radiation 
integral based on body currents excited on the body due to any form of 


5-6 



voltage vector excitation. For scattering problems this excitation is k**. The 
received (image) and transmitted polarizations correspond to a, p = 0 or <p; 

W(k“**) is a standard FFT window weight function such as Hamming, 
Hanning, etc.. 

The downrange k** = k, vector component is centered on the free space value of k. 


. 2w . . 
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*■0 


(4-7) 


Once we have computed the scattered field E“ ,p (k inc ,k rt ) in bistatic k space, equation 5-2 
is the Fourier transform of the k space scattered field. The key is to compute the scattered 
field in terms of It"* vector components that naturally correspond to transformed image 
spatial coordinates. 

The bistatic k space scattered field is computed from the standard radiation integral once 
body currents have been obtained for any given excitation: 


£«•’(£*,£**) - / { A* Pit*) } e*~* dS (4-8) 


where n“ “ # " * is the receiver polarization unit vector, J* " • or * are the body currents 
resulting from plane wave excitation at the angle corresponding to k” 6 and polarization p. 

Approach: The steps required to obtain an image are: 1) Compute body currents due to 
desired excitation, e.g., plane wave for scattering or port voltages for antenna radiation; 2) 
Specify image extent centered over the body, 3) Specify the approximate resolution desired; 
4) Compute the scattered electric field in k*** = (k <k ” ra , k 00- ) space; 5) Pad the data array 
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with zeros to obtain an array (linear or square) that is a power of 2 in length and perform 
the FFT; and 6) Resort the data and display the one, two, or three dimensional image. 

Image extent, resolution, k space bandwidth, and number of steps in a Fourier transform are 
interrelated. Image resolution and bandwidth are inversely related such that fine resolution 
requires large bandwidth: 


A r = 


2 * 
A* 


(4-9) 


Window extent (centered over the body) is an integral number of resolution cells M: 

Window extent * M AR (4-10) 


M is also the order of the Fourier transform, not necessarily a power of 2. The user then 
arbitrarily establishes desired independent parameters and computes the others from the 
above two equations. 

Two factors limit resolution. The first is the number of basis functions per wavelength used 
to represent the current distribution. If too great a resolution is requested, then one images 
the individual basis functions. This is not physical. The second consideration is more 
fundamental. Since the spatial current behavior is inherently on a wavelength scale, it makes 
little sense to specify resolutions greater than X/2. 

In MOM3D the user specifies an exact window extent and an approximate resolution to 
compute the number of increments M in k*“ space. The actual resolution is then 
recomputed using M. The bistatic electric field produced from the currents excited on the 
body is computed at these M k space locations in down range or at M by M k space 
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locations for 2-D down/cross range images. The Fourier transform is then done using a Fast 
Fourier Transform algorithm. Since the FFT requires data arrays that are powers of two, 
the data are padded with zeros to the next largest 2 N . Zero padding smooths the resulting 
transformed image data, however, the original resolution remains unchanged. 

The radiation integral used to evaluate E* requires that Id** be input in terms of rectangular 
coordinates. The fundamental variation for Id'* is (k„ k* k^). We must transform this 
representation back into a rectangular form. If we are computing an image from angular 
location (0,<p), then the rectangular form for k is just the usual vector sum of components. 
Two cases occur. The first is for an "azimuth" or <p directed cross range direction where 
(lq,k # ) are transformed: 


k* = k r sin& cos<|> - sin<j> 

kj * k,sin6 sin<t> + A^cos<t> (4-11) 

k* = k r COS0 


while for an "elevation" or 0 directed cross range direction the transform of (Iq, k,) becomes 


£/ = k r sind cos4» - Jfcg cos<J> cos0 

kf * Jfc r sin0 sin4> + ifcg sin<(> cos0 (4-12) 

k£ * Jt r cos0 - itgSinO 


Lastly, the form for the windowing weight function W is unity at the array center and tapers 
to zero at the array ends. A radial hat function may be used, [14], or a dual weight function 
with a separate function for rows and columns may be used. 
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Cross Range Direction Rotation: For two dimensional images one may choose the cross 
range direction in either the or 0 directions. The $ cross range direction corresponds to 
an azimuth rotation about the Z axis. This would be the typical case for an aircraft 
geometry coordinate system in which the Z axis is vertical. The 0 cross range direction 
corresponds to a polar angle (elevation) rotation about an axis in the X-Y plane that is 
perpendicular to the specified azimuth direction <j>. 

Rotation Center: The image is computed by varying the down and cross range values of k”*. 
Down range direction implies that the corresponding radial direction points toward the origin 
of the coordinate system. This origin is then the center of the image window extent for the 
resulting image. Often, however, the coordinate system for the body geometry has an origin 
that is not necessarily where one wishes the center of the window extent For example, most 
aircraft coordinate systems define an origin at or near the vehicle nose. For imaging, one 
would like to place the origin midway down the body so that the resulting image fills the 
window extent To accomplish this translation shift of coordinates, a user could change all 
the vehicle coordinates in the geometry model file. A simpler approach is to introduce a 
phase shift that corresponds to the desired coordinate translation. 

The k space scattered field for the original or "old" bistatic scattered field radiation integral 
is : 




(4-13) 


When a new reference origin, R^ is defined relative to the previous body origin, then 
the new image transform becomes: 
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£"(£') - I /(£') ' <*•£-> <£ 

O^) - /(f) ««*** dS < 4 - ,4 > 

O - ‘ O 

Therefore a coordinate shift is implemented as a phase shift at each k space scattered field 
value prior to Fourier transforming the data. We do need to remember that the coordinate 
rotation matrix then applies for a given viewing angle: 


‘ 

* 

cos(0) 

-sin(0) 

X 

(4-15) 

yrc, . 


sin(0) 

cos(0) 

y. 



Example results for down and down / cross range are shown in the example results section. 
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SECTION 6 


NEAR FIELD COMPUTATIONS 

The theoiy for computing near zone electric fields is presented. Near fields are computed 
from body current distributions that result from a specified voltage vector excitation and 
system impedance matrix. The theory applies for electric current sources and therefore 
applies for boundary conditions corresponding to a perfect electric conductor or to resistive 
surfaces. When compared to far zone fields, the near zone field computation requires the 
addition of the scalar potential term to the formulation. 

Near surface electric fields are of interest to understand the scattering process and to 
evaluate E field magnitude, phase, and direction on or near a surface. Quantities involved 
are the total field, the incident field, and the scattered field. The computation of surface 
E fields is the first step in computing surface power flow. 

Far field radiation is usually of interest and is computed by most all algorithms (including 
MOM3D). The far field is directly proportional to the transverse components of the vector 
potential and is used to compute radar cross section or antenna gain patterns. 

The near field, in contrast, involves both the scalar and vector potential (conservative and 
solenoidal sources). While the near field expression is perfectly valid for the far field, one 
seldom goes to the trouble to compute the more complicated expression if only the far field 
is of interest 

Far field scattered radiation field lines are solenoidal, i.e., they close back on themselves 
without ending on a charge source while near fields have both solenoidal and conservative 
vector components. On a perfect electric conductor (PEC) surface the field lines originate 
on the induced surface charges. Near a PEC surface the field is dominated by the scalar 
rather than the vector potential and the field can be represented as a quasi static. As the 
field point moves away from the surface, the vector potential term increases its contribution. 
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In the far field the radial field from the scalar potential identically subtracts from the radial 
vector potential component leaving only the transverse components of the vector potential 
as the final contributor to the far field radiation. 

Far field radiation usually involves the scattered field, i.e., that radiated by the surface 
currents. Near field radiation, in contrast, must consider the total, incident, and scattered 
fields, E* = E’ + E\ 

Boundary conditions for a perfect electric conductor (PEC) require zero tangential total 
electric field, i.e., a PEC surface is a short circuit that does not support a tangential electric 
field. Only a perpendicular field can exist at the surface of a PEC E‘ is composed of the 
incident field E' and the field scattered by the body, E*. The physical process by which the 
boundary condition is satisfied is by induced surface charge and current The scattered field 
caused by these sources, E f , is equal in magnitude but opposite in direction to the incident 
field E' on the body surface (tangential component). The total tangential field is zero, E*^ 
= [ E* + E* ]„, = 0 and thus E^ = - E^. 

Scattered Electric Fields: The fundamental theory for near zone scattered field is reviewed 
and the expressions for incorporation into MOM3D are derived. A discussion of the various 
field types and their physical meaning is then discussed. 

Currents and charges induced on a scattering body radiate both electric and magnetic fields. 
A far field EM wave is composed of transverse components of E and H whose magnitudes 
are related by the impedance of free space, Z = E/ H = [/i 0 /€ 0 ] 1 ' 1 - 377 - 120 iz Ohms. 
In the far field, E and H are not independent quantities, we may compute either. 

In the near field, E and H are not related and must be computed separately. In this work 
we compute only the scattered electric field. 

The time and frequency dependence of the fields have been assumed to be time harmonic 
and monochromatic, e** 1 . 
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Electromagnetic theoiy tells us that the scattered electric field is composed of vector and 
scalar potential terms for the solenoidal and conservative components of the vector field. 
When dealing with PEC or RBC boundary conditions, only electric current and charges are 
sources of the scattered field. In the more generalized case where magnetic currents and 
charges are required, additional terms are needed. 

The expression for the scattered field at spatial position r f is: 


£*"(/>) - -;W - V,* 


(6-1) 


where j is the square root of -1, <■> is the radian frequency 2*f, A is the vector potential and 
4 is the scalar potential with the gradient taken with respect to the field coordinates. 

The vector potential A is an integral over the vector surface currents: 


A * l » 0 ffgdS 


(6-2) 


where /x is the permeability and J the surface current. The Green's function is 


g(? r r t ) 


e -fi<r r V 

4 * 1^1 


(6-3) 


where k points in the direction of the scattered field and has magnitude of 2ic/A, and the 
vector R is from the source coordinate r $ to the field coordinate r r , i.e., R points away from 
the source point. 
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The scalar potential has electric charge as its source: 


*(f,) 


— J ogdS 
^0 


(6-4) 


where a is the electric surface charge density. 

Surface current J and charge density a are not independent quantities. They are related by 
the requirement of charge continuity as expressed by: 


V, -J = -j<oa 


(6-5) 


where the surface divergence (derivative) is with respect to source coordinates. Thus we can 
rewrite the expression for the scattered field entirely in terms of the surface current J: 


£“*(/>) « -jkx\ f [Jg + (V,-7)V^/Jfc 2 ]dS (6-6) 


where we have used kq = tofi and toe — k/q where n (=377) is the impedance of free space. 
V, is a derivative operator on the field coordinates and slides through to g(r t ,r s ). The last 
step is to use the relationship that V, = -V, and that the gradient of the Green's function 
on the source coordinates is expressed as: 


W /.0 - -<1 * W ) s 


(6-7) 
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where the vector R = r, - r, is from the source to field point 


Finally, we combine the above to arrive at the expression for the scattered field: 


- -jkti / [7 - (V,-7)(l *jkR)(f r T,)Hkg.f ]gds «-*> 


This is our desired result which is a function of the surface current J. It is evident that this 
expression involves additional computation when compared to that required for only the far 
field. Once J is computed for a specific frequency and excitation the near zone electric field 
is computed with this expression. 

Scattered Field Expression for MOM3D: All that remains now is to express the current 
vector J in the basis function representation used in the surface patch model of MOM3D. 
The triangular couple sub-domain basis function expansion is: 

7 - £ ;,/,(?) <«•»> 

i-1 


where j s is the complex current coefficient determined from the solution and f(p) is the 
vector basis function: 





; S* • +1 ; for ? € A* 


; S" * -1 ; for ? € A" 


( 6 - 10 ) 


6-5 



where p is the vector from the triangle vertex to any point in the triangle, 4 , is the length 
of the qth edge, and A, is the area of the qth triangle. This basis function has the property 
that its magnitude is unity any time p lies on the edge (. The surface divergence of the 
basis function has the form: 



( 6 - 11 ) 


We note that the vector p points from the vertex to the opposite edge of the basis triangle. 
Current flows from the + triangle of the couple to the - triangle (by definition). MOM3D 
geometry defines p as pointing away from the vertex, despite whether the triangle is + or 
Thus we introduce a unity sign term that is either + or - depending on which triangle 
of the couple is being considered (or zero if the triangle edge does not have a common 
neighbor). 

Combining the above expression for J and V*J into the expression for scattered field, and 
replacing the integral with a sum over couples, we obtain the expression used in MOM3D 
to compute the scattered near field: 


r~V,) - -y*n £ E / 1 •£ - (1 <» < M2 > 


f-i *-1 


where the first summation is over all the triangles and the second summation is over each 
vertex. This sum is over all the current couples. 
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The integral over the triangle is approximated using the centroid approximation (similar to 
MOM3D matrix elements), i.e., the integrand is approximated by its value at the triangle 
centroid: 


i / F(n dS » F(f t ) 


(6-13) 


where r e is the centroid. Using centroid evaluation for the integrand yields the desired 
expression for the scattered field: 


= -jk ti Y, £ s 

i-i »-i 



(1 *JkR e ) 

(k* e f ‘ 


8(W 


(6-14) 


where R,. = r, - r w and r v is the centroid of the source triangle. 


Every Thing You Wanted to Know About Complex Vectors But Were Afraid to Ask: The 
vector quantities have components that are complex numbers. Each spatial vector 
component has its own magnitude and phase. Each vector is understood to be multiplied 
by the time dependence factor e***. 

Real world physical quantities, however, are not complex. Complex notation is the formal 
process by which the time delay due to the finite speed of propagation is expressed by phase 
(delay). The physical quantity is always understood to be the REAL part of the complex 
number. For an arbitrary complex vector E with e i **‘ time dependence, the physical value is 
the REAL part: 
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(6-15) 


<•1 

where |Ej| and <p { are the magnitude and phase of each vector component of E and the sum 
is over the three vector components with unit direction e,. The physical vector quantity can 
then be written in rectangular coordinates as: 

a (6-16) 

E fkyncoi * lEjcosOk + wr)* + |E y |cos(<|> y + or)y + lEJ cos(<j> t + wf)z 


where we have used the identity e* = cos(x) + j sin(x) and where (xj,z) are the standard 
unit vectors for an orthogonal representation. 

We see that this complex time harmonic vector changes direction in space as a function of 
time. We could display this time dependent behavior by adding «t time increments between 
[0,2it] to create a dynamic moving display of the solution. 

However we also can investigate the two principal directions of the vector by taking the 
REAL or IMAG parts of the complex vector. To illustrate, we first set time to zero, ot = 
0, obtaining: 

* |£jcos(4>, + 0)je + \E y \cosW y + 0)f * |£ z |cos(<fr z + 0)£ 

(6-17) 

- *E X * + *E y f + 8t E t t 
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while if we set time to - 90 degrees, &>t = -ic/2, we obtain: 

9 | 90) i* |£ y |cos(<|>,-90)j> ♦ |EJ cos(<|> z -90)£ 

« \E X | sin(<}> x )i + |£,|sm(4> y )5> + \E t \ sinMJt 
= 3E X * + 3E y ? + S£ z £ 

* 3^0 

(6-18) 


The REAL and I MAG components of the vector E represent the two basis vectors 
(magnitude and direction) that describe the time varying vector. The IMAG part lags the 
REAL part by 90 degrees. A plot of REAL(E) and IMAG(E) will show two snapshots in 
time of how the time varying vector behaves. Every other point in time will be linear 
combination of these two vectors. 

We now review the idea of time average quantities. Sensors seldom measure actual time 
variation, rather a time average value is "sensed." The scalar root mean square (RMS) time 
average of a vector quantity is obtained by averaging the real components: 

E^ * < 8t(£ e j9 ') ‘X(£ eJ*') > m ( 6 * 19 ) 


where < > indicates an average (integral) over o>t = 0 to 2u divided by 2k and we have 
taken the dot product of the two REAL vectors. Each component of the dot product has 
a term that is proportional to: 
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Ifip 

2 * 


( 6 - 20 ) 


/ cos*(4>| + «f) 

where * indicates the complex conjugate for each vector component i. Thus the RMS time 
average of a complex vector is simply obtained as: 



( &_ »r V* 

V5 


(6-21) 


where we recognize that v ^/2 = 0.707 is the root mean square (rms) average of a sinusoidal 
varying quantity that has a unity peak value. 


While the above discussion was for arbitrary complex vectors, we now review the specific 
vectors required in this work. They are the incident field E', the scattered field E 1 , and the 
total field E*. 

The incident field is specified by a plane wave with polarization, wavelength, and direction 
of propagation: 




( 6 - 22 ) 


where u is the polarization unit vector, k is the direction of propagation, and Eq is the scalar 
amplitude. In MOM3D we take Eq = 377 which represents a unity value for the associated 
magnetic field This is done so that the currents are then normalized to Hq for 
comparison to Physical Optics values, 2nxH'. 
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The scattered field E* is that due to currents on the scattering surfaces given by the 
previously derived expressions. 

The total field E' is the vector sum of the incident and scattered fields, E* = E 1 + E\ 

Near the scattering body, we mostly are interested in the total field since that is what a 
measurement probe (either dynamic or average) would measure. In the far field, a 
measurement probe, e.g^ a radar receiver, would still measure the total field, but since the 
transmitted field has been turned off (and is in the opposite direction) when the scattered 
energy arrives back from the target, the measured quantity is just die scattered field E*. 

Formulation When the Field Point is Near Surface: When the field point for computing the 
scattered field is within one triangle dimension of the surface, the centroid approximation 
breaks down due to the 1/R and 1/R 2 singularities in the formulation. The physical field is 
of course finite being due chiefly to the surface perpendicular field from the surface charge 
density a. This section discusses the analytical approach used for computing the scattered 
field when the field point is very near the surface. 

Several approaches could be developed. One could better approximate the singular integral 
with either a closed form analytical integration (if it could be derived!) or develop a surface 
numerical integration scheme. But, since only those few field points very near the body 
surface are affected by the singularity, it is not worth while to develop elaborate schemes to 
"fix" the problem. Thus the following approximation is used. 

When the field point is on the surface with R * 0, the singular self patch field is its 
perpendicular value. This same vector value is then used to approximate the near surface 
self field as the field point moves away from the surface to a distance no greater than a 
triangle dimension. 

Electromagnetic theory and boundary conditions tell us that the surface parallel and 
perpendicular fields are (for a general surface): 
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(6-23) 

/?!*’ ^ ( Physical Optics Assumption ) 

it 


where ( o, J ) and ( o', M ) are the local surface electric and magnetic charge density and 
current For our MOM3D analysis, we have only resistive and perfect electric conductor 
boundary conditions, therefore equivalent magnetic sources do not exist thus: 


By* at £y* = 0 

2«o 1 


(6-24) 


We thus make the assumption that very near the surface, the fields are the same as on the 
surface. We will approximate the self field with the above expression whenever R = R, - r, 
is less than approximately one linear triangle dimension and use the centroid approach 
whenever R > d. 

The self patch parallel field is zero. The perpendicular component is developed by 
expressing the surface charge density in terms of the surface divergence of current: 
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&•* , A , a 

A 2€ 0 2k 

where V*/ » -;'uo , and ue /i) 

flimensioMl Cfedt : 2*2. SSE*!” . 

1/m m meter 


In MOM3D, the triangle couple basis functions result in the following expression for the 
surface divergence: 


V-/ = +-?L or 

a; a; 


(6-26) 


The resulting expression for the i th self patch perpendicular field is thus: 


gacatjtff _ 


M . 

± 2k A* Ji 


(6-27) 


This egression is used for evaluating near surface patch fields where the 1/R singularities 
dominate. For field points further away, the previous approach using the centroid 
approximation to the integral is used. 
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SECTION 7 


MONOSTATIC and BISTATIC RADAR CROSS SECTION 

Monostatic (backscatter) and bistatic radar cross section may be computed in MOM3D. The 
body is excited via a plane wave from a specified direction in space, k**, and the resulting 
body currents are computed. The far field radiation resulting from these currents is then 
computed. For backscatter RCS, the scattered field is computed at the same angular 
location as the source radiation. For bistatic RCS the scattered field is computed at angle 
cuts specified by the user. The required theory for computing radar cross section is 
presented. 

Polarization Scattering Matrix: A complete electromagnetic scattering description is 

contained in the polarization scattering matrix that relates the two possible transverse 
incident field vectors to the two possible scattered field vectors. Only two components of 
E are allowed since Maxwell’s equations require that the E field be transverse to the 
direction of propagation. The polarization scattering matrix relates the incident field vector 
to the scattered field vector: 


B"* 


(7-1) 


where S is a 2 by 2 dyadic (tensor) that relates the components of E”* to E": 




s ®* 6 s 9A 


£®**® 



s*fi 




(7-2) 
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The diagonal terms are the two independent co-polarized scattering elements while the off 
diagonal terms are the two cross polarized scattering terms. The scattering matrix 
completely represents the scattering properties of a target The four terms of S are complex, 
each with amplitude and phase. The elements of S are related to RCS: 


^ 4 nr 2 


(7-3) 


E“ and E“* are independent functions of angle, thus the elements of S become functions 
of the source and receiver angular coordinates, 




(7-4) 


All possible electromagnetic scattering information is contained in the polarization scattering 
matrix. Huynen [15] and others have attempted to characterize (and therefore identify) 
target properties such as size, orientation, symmetry, de-polarization, and characteristic 
angles, etc. from the information contained in S. 

MOM3D has the potential to compute the complete scattering matrix for both linear and 
circular polarization. 

Radar Cross Section Definition: We generally do not work directly with the polarization 
scattering matrix S, we deal with the components of S as expressed by the co- and cross- 
polarized radar cross section. Radar Cross Section (RCS) is a measure of the power 
scattered by a target normalized to the power density incident on a target. Except for 
incident and received polarization, RCS is defined to be independent of specific radar system 
parameters such as transmitter power, receiver sensitivity, range to target, etc.. 
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The definition of RCS assumes an incident plane wave on the target (which physically does 
not exist!). The scattered field is then computed or measured far from the target to avoid 
near field effects. Radar cross section is formally defined as the far field scattered power 
density normalized to the incident power density at the target: 


,«.P 


Lim 4icr 2 — ; P « E 2 or H 2 

2 m 9 




r - « 
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where a and p represent the polarization (6 or <p) of the scattered and incident radiation. 
RCS, while independent of specific radar system parameters, is a function of: target 
geometry; frequency ( or wavelength ); incident polarization; received polarization; and 
the angular position of the source and receiver: 




(7-6) 


The power density for an electromagnetic wave is the sum of the E and H field components. 
For plane waves the E and H fields are related via wave impedance. Thus wave power 
density can be expressed entirely in terms of either E or H. Once the body currents are 
obtained, RCS may be computed via either E or H since in the far field they are related by 
q = E / H. Typically we compute the scattered E field assuming E to is unity. However, is 
it of more interest to have H"* = 1 so that we may compare resulting currents to normalized 
physical optics values, J = 2n x H“. 

Radar cross section computation requires a scattered field, which in turn requires a body 
current distribution as excited by an incident plane wave. The formal theory for the plane 
wave voltage vector was presented in Section 2 and 3. Once a body impedance matrix has 
been computed and LU decomposed (or the inverse found), body currents corresponding 
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to the voltage vector are determined. The field radiated by these currents is computed using 
the row measurement vector (Sections 2 & 3). 

The voltage vector used to compute the body currents is a function of the angular 
coordinates and polarization of the illumination source while the row measurement vector 
is a function of the angular coordinates and polarization of the receiver. 


R‘ - £•(&**, <|>**) ; K p - F p (0‘* c ,4> toe ) 


(7-7) 


Radar cross section can then be defined in terms of the row measurement vector and body 
currents: 


■* = Ifv'l 1 

4rc 


(7-8) 


where we have assumed E” c = r| = 377. For the co-polarized case, a = p = 0 or 4>. The 

/ 

current vector j is formally expressed as the inverse of the system impedance matrix times 
the excitation voltage vector: 


p - [zr 1 p p 


(7-9) 


Monostatic Case: Backscatter RCS is computed by evaluating the row measurement vector 
at the angular location of the plane wave source illumination. Our Galerkin weight function 
approach yields the same computation form for the voltage and row measurement vectors. 
Thus the backscatter row vector is equal to the voltage excitation vector used to compute 
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the body currents, R , (0“ a ‘,^**) = V , (8“ c ,<fr” e ). Therefore the backscatter RCS subroutine 
does not need to compute the R row vector for the co-polarized case. 

Bistatic Case: Bistatic RCS is the case where one is interested in the scattered radiation at 
angular locations that are not the same as the illumination source. Bistatic RCS requires 
specification for both illumination and receiver angles. Two cases arise: 

1) The illumination angle is fixed relative to the body while die receiver location is 
moved (angle "cut”); and 

2) The illumination source and receiver are fixed relative to each other while the 
body is rotated. 

In case 1) only one set of currents is computed for the body and die bistatic RCS is 
computed by re-evaluating the row measurement vector for each receiver angle of interest 
This is the most common "analytical" case and is implemented in MOM3D. 

In case 2) the body currents must be computed for each new illumination angle. For this 
case, a new voltage and row measurement vector must be computed for each angle and V 
is no longer equal to R. This is the most common "experimental" case and corresponds to 
rotating a target with fixed illumination and receiver angular locations (with constant angle 
between). This case is not presently implemented. 

In each case the row measurement vector must be computed at the new receiver angle. This 
is done in MOM3D using the subroutine that computes the voltage vector V. 
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SECTION 8 


CURRENT COMPUTATION 

In efforts to understand the radiation and scattering mechanisms it is often desirable to 
display the current distribution on the body as caused by either in incident plane wave (as 
in a scattering case) or by local surface "port" excitation (as in an antenna problem). This 
section discusses a few of the issues required to compute the currents, spatial and temporal 
average values, and the two independent time values. Current "modes" of interest might be 
those corresponding to: physical optics; surface traveling waves; creeping waves; edge waves; 
leading edge diffraction; trailing edge diffraction; and multiple bounce. 

Once the body impedance matrix is computed and solved and the voltage vector specified 
(either as an incident plane wave or as "port" excitation), the unknown current vector j is 
formally computed as: 


P * [Z] ' V* 


(8-1) 


The complex vector body currents are then given by the surface basis function expansion: 
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Couple Spatial Averaged Current: The basis functions are functions of position p within 
each triangle. The triangle hat basis functions can be averaged over the couple (two 
triangles) as suggested by the results of Section 2. On a couple basis, the current average 
is given by: 
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where ri* is the centroid location of the + and - triangles of each couple. This is the 
quantity that should be considered as representing the surface currents. The location of this 
averaged vector can be taken to be either end of the vector, ie., triangle centroids, or the 
midpoint location, 0.5 * ( p e+ + p c ‘ ). If the two triangles making up the couple are not 
co-planar, the midpoint is not in the plane of either triangle. 

Triangle Currents: The number of 

triangles making up the geometry is not 
equal to the number of current couples. 

Often for ease of display it is convenient to 
show currents at the centroid of each 
triangle. This can be done simply by 
performing a vector sum of the current 
couples associated with each triangle, 

Figure 8-1. This sum is over one, two, or Figure 8-1 Triangle current can be 
three couples depending the number of considered as the sum of couple currents 

common neighbor edges. Edge triangles 

have one free side ( two couples ) while a tip triangle has two free sides ( 1 couple ). The 
triangle spatial average complex vector current is thus: 
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where J q "’ e is computed from the previous expression. 
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Time Varying Vectors and Surface Current Animation: Surface currents are time varying 
vectors. The physical quantity is obtained by taking the real part after multiplying by the 
time phasor e***. At each couple location the current vector sweeps out a trajectory on the 
sur f ace as o>t varies between 0 and 2a. Two principal vector directions (basis vectors) from 
which all other values can be represented are the vectors at «t = 0 and ot * -a/2. These 
two vectors are the REAL and IMAG parts of J: 
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(8-5) 

/~((or--a/2) - 


The surface current vector as a function of time is a linear combination of the real and imag 
basis vectors: 
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A dynamic graphical display of currents would simply vary wt over [0,2a]. One could display 
die dynamic vectors with magnitude shown as relative vector length or color. Or one could 
display the magnitude using a color scale without showing the vector direction. These 
displays would represent the dynamic currents induced on the body due either to an incident 
plane wave ( as for RCS problems ) or to local surface port excitation (as for antenna 
problems ). For plane wave excitation one also could display the incident wave (forcing 
function) as it passes over the geometry. 

Time Average Currents: Besides a spatial average, there is also a time average to consider. 
The time average is a stationary scalar quantity (not a vector) obtained by averaging the 
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REAL part of { J e* - * } over the interval [0,2k]. The root mean square (RMS) time average 
is computed in the same manner as discussed in Section 6 for near fields: 
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This stationary scalar quantity could be graphically displayed using a color magnitude scale 
for either the couple currents or the triangle based currents. 
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SECTION 9 


ANTENNA FORMULATION 

MOM3D can be used to compute antenna gain patterns. Since the impedance matrix 
represents the electrical interaction of the body with itself independent of the excitation 
voltage vector, we can use a voltage vector corresponding to a local surface "port" excitation 
to excite body currents. Radiation from these currents is the antenna pattern. This section 
will discuss the voltage vector, antenna gain, radiated power density, power input, and input 
impedance. 

Geometries operated as antennas have voltage forcing functions that are local. These can 
take the form of a coax feed from a transmitter or simply an open ended waveguide 
terminated on the surface over which there is an E field. The antenna problem does not 
depend on the nature of the transmitter, it depends on how the structure is excited: voltage 
magnitude, phase, and vector direction. This excitation induces currents on the body 
according to Maxwell’s equations and boundary conditions. 

Voltage Vector: For the antenna problem, the local surface voltage excitation occurs at any 
one of the surface couple basis functions. Each excited location is called a "port." For most 
antenna problems only one port is excited. But there may be multiple excitation ports on 
a body, each with its own magnitude, phase, and direction^ The voltage vector forcing 
function takes a simple form: 
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Each excited port element of V is computed as before with W = f: 
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where we have assumed the applied port electric field E* is constant over the couple (two 
triangles), the expression for f** has been used, and Ar c * is the vector from the + to - 
triangle centroid. In MOM3D we take E to be in the same direction as Ar so that the 
applied voltage across the couple V tenB is: 


v urm * f 
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In MOM3D we input the value of V, directly so that the terminal voltage is expressed as: 
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Once Vj has been input for one or more excited ports, the body currents are computed in 
the normal manner: 
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The far field radiation at specified angular direction and polarization a = 0 or 4> is 
computed in the usual fashion in terms of the row measurement vector for the desired 
polarization: 
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Antenna Gain: The antenna pattern or gain is computed in terms of the far field radiation 
power density normalized to the input power to the antenna (assuming the input power is 
radiated isotropiclly and that we have unity efficiency): 
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where has units of watts while P™ 1 has units of power density, watts/meter 2 , i.e., is the 
Poynting vector of the far field radiation. The astute reader will notice the similarity of this 
definition to that for radar cross section. 


In the far field where E / H = ti = 377, the power density is the REAL part of the Poynting 
vector 


P 


rad 


0.5 8t(fx#’) 


I |2 

2t, 


(9-8) 


where * represents the complex conjugate and the factor of 1/2 results from using the root 
mean square value for the field quantities. 
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The power input to the antenna is the sum of input powers (real part) of all the excited 
ports: 
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The current at each port is the computed current density j (amps/meter) times the common 
edge length: 
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The terminal voltage was previously given in terms of the port voltage vector so that the 
input power to the antenna is computed in terms of the input voltage port value and 
resulting current coefficient j s : 
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The final expression for antenna gain relative to isotropic is then 

G«(0,$) * ^3. 1^ /l 2 (9-12) 
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This is the antenna pattern for polarization a = 8 or q>. When expressed in decibels it is 
known as dBi, i.e., dB relative to isotropic. A half wave dipole has a figure eight gain 
pattern with broadside maximum of approximately 2.1 dBi. 
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Input Impedance and Admittance: Another quantity that is often of interest is the input 
impedance of the antenna at each excited port. This is obtained directly from the definition 
as the ratio of port voltage to current: 





(9-13) 


Z has both amplitude and phase. The real part is the radiation resistance while the 
imaginary part represents energy storage, either capacitive or inductive (current lagging or 
leading). A half wave dipole at resonance has an approximate input impedance of 72 ohms. 
At resonance, Z is often completely real, i.e., zero degrees phase. 

The input admittance is just the reciprocal of input impedance: 



(9-14) 
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SECTION 10 


RESISTIVE BOUNDARY CONDITIONS 

Resistive Boundary Conditions (RBC) have been implemented in MOM3D. This section 
describes the ohms per square concept for thin resistive materials and its extension to the 
more general case of thin dielectric material. Resistive Boundary Condition (RBC) is then 
presented followed by the more involved theory for Impedance Boundary Conditions (IBC) 
and the approximate IBC case. This is shown to be the same as the RBC case, therefore 
the present RBC MOM3D can be used where the approximate IBC hold. 

Resistive Sheet Concepts: Resistive Boundary Conditions (RBC) are typically used to model 
thin sheets whose impedance does not change with angle of incidence. Such sheets typically 
are thin resistive layers, usually frequency independent, or thin dielectric layers (which are 
frequency dependent). Thin layers are characterized by a complex sheet impedance known 
as Ohms per Square (OPS) that has magnitude and phase, i.e., characterized by energy 
dissipation and reactive storage. 

Ohms per square characterization of thin resistive layers is derived directly from bulk 
material resistivity characterizations when the thickness x is allowed to become a thin layer. 
The resistance of many materials is characterized by resistivity p ohm-meters or by its 
inverse, conductivity a = 1/p, whose SI units are Siemens per meter. S/m, or mhos/m (mhos 
= 1/ohms). 

The OPS concept is derived by considering a rectangular block of resistive material. The 

resistance between parallel faces of the block, R increases with spacing f, and decreases 

with area A, Figure 10-1,: 
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where t is material thickness and w is 


width. 

When thickness t becomes small, as for a 
thin sheet, then the resistance between 
terminals is expressed in terms of 
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Figure 10-1 Ohms per square concept 
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where has the definition 
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Rop, then becomes the electrical characterization of the material rather than specifying 
conductivity o and thickness x. The resistance R^ between sheet ends becomes. Figure 
10 - 2 , 
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If one measured R tCTB between the ends of 
a square sheet, where f = w, then R,^. = 
Rep,, therefore the name Ohms per Square. 
While we generally are interested in the 
values of R^ at microwave frequencies, 
experience has shown that conductivity is 
usually frequency insensitive, therefore a 
DC measurement of R^ is representative 
even at microwave frequencies. Typical 
values range from 0 (conducting surface) to 
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R term - R(ops) */w 
Figure 10-2 Ohms per square (OPS) 


several thousand ohms per square. A unique value is q/2 = 377/2 = 188.5 ops where 
Physical Optics currents are 1/2 that on a PEC surface. 


Dielectric Thin Sheet Impedance Concepts: The generalization of thin resistive sheets, i.e., 
sheets characterized only by dissipation characteristics can be extended to thin sheets that 
also may have energy storage characteristics, such as thin dielectric material, e.g., fiberglass 
sheets. The concept of thin means that the electrical thickness is small compared to a free 
space wavelength: 
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When this holds, the dielectric polarization currents perpendicular to the sheet direction are 
very small compared to the "in plane" currents. In plane energy storage and dissipation can 
then be characterized as ohms per square that is complex, i.e., amplitude and phase. 

A general thin sheet expression is given by [16]: 


10-3 




(10-6) 


^o/o * 


SUL 


k(e r -l)x 


-760 A, 
(€ r -l)T 


where t) is the impedance of free space, (jijt 0 ) l ' J * 120* * 377 ohms, k is the wave number 
2x/A., k is the wavelength, and €, = €'- je' is the relative dielectric constant The loss part 
of c, is c* and when conductivity is the principal loss mechanism, which is the case for most 
materials of interest, it is 
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where <•> is the radian frequency, a the material conductivity, and €q is the permittivity of free 
space. When e" dominates the relative dielectric constant, as it does for resistive materials, 
then the thin dielectric sheet impedance reduces to the resistive sheet expression (using (i>€ 0 

= k/n), 
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Thin dielectric layers, at high frequencies, can have low impedances. This must be 
considered when backing thin R card material with dielectric support surfaces such as 
fiberglass ( - 600 ops at 10 GHz for e r =4, and x = 0.040'). 

Resistive Boundary Conditions: In the original MOM3D code, perfect electric conductor 
(PEC) boundary conditions were imposed, i.e., the total tangential electric field was set to 
zero. 
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where E" is the incident electric field and E* is the field scattered by the body. E scattered 
is an integral over the surface currents and is written in operator form as E* = - L(J). The 
PEC boundary conditions then take the form 
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A solution is obtained by taking the inner product with a set of weight functions and 
expanding the current density J with a set of sub-wavelength basis functions. A matrix 
equation results for unknown current coefficients with the incident wave voltage vector being 
the solution fordng function. 

Resistive Boundary Conditions for thin materials follows Ohms law, namely that the total 
tangential electric field is equal to the current density J times the resistance R, 
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RBC boundary conditions for the EFIE take the form 
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As with the PEC approach, a solution is obtained by taking the inner product of each side 
with a surface vector weight function W, 




(10-13) 


We see that the RBC differs from the PEC only by the term = < W, R J > that is 
added to the diagonal elements of the PEC body impedance matrix. 

The explicit expression for is obtained by using the triangle couple expansion functions 
defined in Section 3. The Galerkin weight functions require W = f. We then use f "* to 
compute Z^: 
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where R^ is the average of the value assigned to the + and - triangles of the couple. 

Impedance Boundary Condition Theory: This section outlines the more general Impedance 
Boundary Condition (IBC) and its approximate form that reduces to the RBC. 

Interest in the IBC is for application to perfectly conducting surfaces coated with a thin layer 
of magnetic (mag) radar absorbing material (RAM). Mag RAM differs from thin resistive 
layers in that performance is due to phasor cancellation as well as energy dissipation. A mag 
RAM coated plate is a Dallenbach layer, Figure 10-3. 
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IBC apply only for high index of refraction 
materials where the path length inside the 
material is always normal to the surface. 
When this occurs, the path length does not 
change with incidence angle. Figure 10-3. 

In the most general case, Maxwells’ 
equations must be solved for both electric 
and magnetic currents, J and M. The 
standard definitions for these currents are 
in terms of tangential surface magnetic field 
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Figure 10-3 Dallenbach layer 


H and electric field E, 


J • + fixff™ 

M = - n x S™ 


(10-15) 


where n is the local surface normal. For PEC surfaces there are no magnetic currents since 
E 70 *^ = 0. For general impedance surfaces both H and E tangential fields can exist and we 
must solve for each. 

Maxwells' equations in the Stratton Chu integral form for closed surfaces express the 
scattered E’ and H* fields in terms of J and M source currents. The EFIE takes the form 

[7]: 

ST - -/ [ycopjg + tfxVg - 0i-£)Vg] ds ( 10 - 16 > 


where g is the Green's function, exp(-jk-R) / ( 4icR ). The corresponding MFIE is [7]: 
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It is customary [17] to express the boundary conditions in terms of L and K integral 
operators: 


i™ - i* - 1(7) + K(M) and 
H™ - r-f(J)-I(AO/ti! 


(10-18) 


where i) 0 = 377. The operators L and K have the definition [17]: 
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where X is either J or M. 

From this general formulation, we see that an arbitrary impedance surface requires solving 
for magnetic and electric currents. 

When one assumes that the surface impedance does not change with viewing angle, i.e., that 
the index of refraction is n » 1 such that phasor cancellation path lengths do not depend 
on look angle, then one can apply the impedance boundary condition. 
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The impedance boundary condition relates the magnetic current M to electric current J via 
the surface impedance Z* 
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which implies that [17] 


M - -AxE T * - -Z s AxJ (10-2D 


When IBC conditions apply, it reduces the unknowns from both J and M to only J. This is a 
very significant reduction in the complexity of the analysis and computer algorithm. 

As indicated in [17], the IBC approach leads to physically correct solutions at interfaces 
where the refractive index of the surface layer is much greater than unity, n = (e/t r )‘ ,J » 
1, and where the surface impedance can be expressed in terms of the surface material c, /*, 
and thickness. 

Surface impedance is constant with viewing angle if and only if the index of refraction is 
much greater than unity, n » l(see Figure 10*3). When this is the case, specular energy 
incident on the layer travels perpendicular to the surface when inside the material layer. 

•v 

Using the IBC to eliminate the magnetic current M, reference [17] obtains: 
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Reference [17] solved this equation for Bodies of Revolution (BOR's) by computing the L() 
and K() operators for the BOR coordinate system and Fourier modal basis functions. 

A similar approach could be done for the MOM3D code coordinate system and basis 
functions if we computed the K operator in addition to the L operator matrix elements 
presently computed. Such a solution is clearly a large undertaking and is beyond the scope 
of the present effort Several concerns, however, are the requirements for closed surfaces 
required for the MFIE and the vector representation of n x J for the triangular patch basis 
vectors. 

Approximate Impedance Boundary Condition Theory; If we require that Z, < < 377 and if 
we approximate the K operator (which is a Fredholm integral equation of the second kind) 
to just its self term, then we can obtain the approximate IBC equation that has the same 
form as the RBC case. The K(X) operator is approximated as: 

= 1 + J JfxVg dS » - (10-23) 

2 tv 2 


This approximation is the tangent plane approximation, valid for surfaces whose radii of 
curvature are large compared to A, or is valid for planar surfaces away from edges. This is 
identical with the Physical Optics approximation where the surface field is due only to local 
sources and does not have contributions from other surface regions. When this 
approximation is made with the requirement that Z, « q 0 ” 377, the approximate IBC 
equation is obtained 
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This is the same form as the RBC. Thus the approximate IBC is the same as the RBC that has 
been implemented in MOM3D. 

Summary: The resistive boundary condition expressions have been derived. The RBC can 
be applied to thin resistive or dielectric sheets whose surface impedance is not a function 
of illumination angle or polarization. 

Impedance boundary conditions were reviewed. The material coating requirement to apply 
the IBC is that the index of refraction in the coating be n >> 1 so that the surface 
impedance is not a function of illumination angle. 

An approximate form of IBC was developed that required that the surface impedance be 
Zjuf « 377. For this case, the approximate IBC formulation is the same as the RBC 
formulation. 
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SECTION 11 


SYSTEM MATRIX SOLUTION 

The motivation for this section is to remind ourselves of matrix solution considerations that 
are a fundamental requirement for method of moments solutions. Topics discussed are the 
system matrix, its dependence on linear body size, factors leading to matrix ill conditioning, 
and a discussion of the matrix condition number. 

System Matrix: The impedance matrix represents the electrical interaction of the body with 
itself at the specified frequency (wavelength) according to Maxwell's equations and boundary 
conditions. All electromagnetic interactions are accounted for: physical optics (specular, 
i.e M angle of incidence = angle of reflection and end region), leading and trailing edge 
diffraction, shadowing of one surface by another, surface traveling, creeping, and edge wave 
phenomenon, multiple bounce(s), etc.. 

The system matrix is a powerful representation of the EM body response since it is 
independent of the nature of body excitation. Once the matrix is computed and solved via LU 
decomposition or inverse, we then obtain surface currents for any type of excitation, e.g., 
RCS scattering or antenna gain patterns. RCS scattering problems can specify polarization 
and angle of illumination for plane waves, or for spherical or even cylindrical waves provided 
one uses the corresponding voltage vector. Or the same matrix can be used to obtain 
antenna patterns with one or more localized regions of the body excited, each with user 
specified amplitude and phase. Once the system matrix is computed and solved, we can save 
the matrix for re-use for any type of EM problem. 

The method of moments approach with excitation independence is contrasted to the 
differential equation approach to solving EM problems such as the Finite Difference Time 
Domain (FDTD) method. The FDTD, while in principle solves for all frequencies at once 
if an impulse excitation is applied, must be re-solved for each excitation. Thus to obtain an 
RCS backscatter plot at 1 degree increments, 360 new and different FDTD solutions must 
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be computed. This is in sharp contrast to the MOM approach where one does not compute 
a new Z but only a new voltage vector and currents. 

The limitation of any MOM code, however, is the memory storage required for the matrix 
elements and the solution time to solve. The number of unknowns for 3-D surface problems 
is proportional to linear body size squared (L 2 - area). Memory requirements vary as the 
fourth power of body size, N 2 « L 4 , while solution time increase as the sixth power of body 
size, N 3 « L 6 . 

We have been able to increase body size as computer resource capabilities have grown. The 
L 6 increase in CPU time implies that for every factor of ten increase in computer speed, the 
corresponding body size capability grows by 47 percent ( 10 1 * = 1.47 ). Similarly, every 
factor of 10 increase in memory increases body size by 78 percent (10 ,/4 = 1.78). Continued 
growth of serial computer architecture will not translate into appreciable growth in body 
size. For present architecture one must wring out all possible efficiencies such has been 
done in MOM3D by using matrix and body symmetries to reduce memory and CPU times. 

Significant increases in body size will most likely require parallel computer architecture with 
vastly improved throughput, or we will need a revised formulation of the problem to reduce 
significantly the number of unknowns. 

MOM3D has treated the actual matrix solution as a black box operation by using the 
LINPAK library of matrix algorithms [9]. Other libraries are available and a motivated user 
should feel free to implement. 

The practicality issues and virtues of LINPAK are: 

* Routines for matrices that are complex; 

* Routines for symmetric matrices with reduction in run time and storage; 


11-2 


* Routines that compute the matrix condition number to check the "quality" of 
the solution; 

* Routines written in FORTRAN transportable code without common blocks 
(data is passed via subroutine arguments); 

* Routines which build on the optimized low level Basic Linear Algebra 
Subroutines (the BLAS) as fundamental building blocks; and 

* LINPAK is in the public domain. 

Matrix Condition Number The concept of matrix inverse implies that the matrix is not 
singular, i.e., the matrix determinant is non zero. Matrix singularity implies that the set of 
equations is not linearly independent. The method of moments matrix formulation arises 
by expanding the unknown currents with a set of N basis functions with unknown 
coefficients. The boundary conditions are then applied to N independent locations on the 
body. The solution of this coupled set of equations is then obtained. Except for internal 
resonances, the matrix is never singular (in theory). While the matrix determinant may never 
be identically zero, it may be numerically small enough to cause unstable solutions. Thus 
it almost becomes mandatory to obtain a quality check of the solution. This can be done 
by choosing LINPAK routines that compute the matrix condition number (which typically 
requires a 10% overhead in computation time). 

Factors which can lead to numerically singular matrices are: 

* Too many samples per wavelength such that the body becomes over sampled 
with each resulting equation becoming less unique (independent). This can 
occur when one tries to sample geometry variations on a finer scale than A./10. 
Electrically small bodies suffer from this effect; 
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* Interior body resonances that occur at specific frequencies. This can be 
overcome by using a combined field formulation but at the expense of more 
complexity with many additional unknowns. Since this effect is always 
narrowed banded in frequency, it makes more sense just to be alert to this 
possibility and not to attempt computation near these few frequency locations. 

* The number of unknowns is large and the precision of the computer 
arithmetic insufficient; 

The matrix condition number k is a measure of the sensitivity of the solution to errors in 
the matrix and/or the right hand side forcing function. The following discussion is adapted 
from the LINPAK Users’ Guide [9] with the notation changed to that of our MOM code 
problem. The general MOM matrix equation and solution is: 
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Errors in the matrix Z due to numerical round off or formulation and/or in voltage vector 
V of magnitude e 
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may lead to possible relative error in the current solution J of 
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The LINPAK routine for obtaining the LU decomposition of the complex symmetric 
matrices is CSPCO. This routine computes a parameter RCOND that is an estimate of the 
reciprocal of the matrix condition number, 1 /k. If RCOND is approximately 10^ then the 
elements of J can usually be expected to have d fewer significant figures of accuracy than 
the elements of Z. If RCOND is so small that in floating point arithmetic it is negligible 
compared to 1.0, then J may have no significant figures. On most computers this condition 
may be tested by the logical expression 

i 

[ (1.0 + RCOND) .EQ. 1.0 ] 

When this expression is true, the matrix is considered to be "singular to working precision." 
As a special case, if exact singularity is detected, RCOND may be set to 0.0, Reference [9]. 

MOM3D computes the condition number as k = COND = 1 / ( RCOND + 10' M ). Low 
values of COND suggests a reasonable non singular solution while a value of 10 +3 ° shows 
a useless singular solution. While we do not have definitive threshold values for poor 
solutions, if COND is large and the resulting current solution looks suspicious, then the 
solution most likely is not correct If the computer has 6 digits of precision, and COND is 
greater than 10* 6 , i.e., RCOND < 10 -6 , then the solution for J may not have any significant 
figures. 

Matrix solution algorithms usually produce a result. It is up to the user to figure out if the 
result is garbage or is a good solution. 
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SECTION 12 


MODELING ISSUES 

The objective of this section is to introduce several modeling related issues inherent when 
using the triangle basis functions introduced by [1], particularly the ZIG-ZAG nature of the 
basis function representation and lack of symmetry. 

The method of moments approach to modeling electromagnetic current behavior requires 
two critical steps. First one must choose basis functions that can adequately represent 
currents and charges on the geometry for the specified excitation. Second, the analyst must 
develop an algorithm that adequately models the body-body interactions using the chosen 
basis function representation. If either of these steps are not adequate, then the 
computation model will not produce physically correct results. The quest for correct results 
however, often times leads to inordinate computation complexity and/or computer resources. 
Thus we need to decide when our modeling process is good enough to answer the questions 
we pose. 

Balance in the modeling processes is required. The basis function choice and sample density 
is inherently an approximation to the surface currents. Since we need both adequate basis 
function representation and an EM prescription, we cannot make up for poor basis function 
choice and density by developing a precise EM prescription. Still, sometimes a higher 
density of basis functions can compensate for a poor EM prescription. 

Surface Current Characteristics: Surface current spatial variation inherently is scaled to the 
free space wave length. It is vector in nature, i.e., has directional characteristics. The spatial 
variation X sets the sampling requirements necessary to model a continuum with a set of 
discrete basis functions. Typically we choose 7 to 10 samples per 2 as sufficient to model 
this spatial change, i.e., 98 to 200 samples per X 2 of surface area (2 basis vectors required 
for surface). But sometimes the geometry or electrical characteristics change faster than 
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the wave length scale and we model this behavior. Current variation near edges, depending 
on the voltage vector forcing function, may have a square root singularity in its behavior. 

Current vector direction depends strongly on excitation, i.e., current flows in the direction 
of the applied electric field. Yet, structure also influences direction. For example, edge 
diffraction behavior usually involves currents either parallel or perpendicular to the edge. 
Surface wave reflection (traveling and edge wave) and how we model the edge or tip can 
influence the current direction and magnitude. Edge wave currents are coupled strongly to 
the edge and "reflect" from the tip. Thus a vector basis function representation that is 
completely general to allow for any type of surface current flow for any possible excitation 
function must allow for many directional possibilities. 


Triangle Basis Function Considerations: The surface triangle couple basis functions 

introduced by [1] and used in MOM3D have certain desirable characteristics. The triangle 
nature allows for modeling doubly curved surfaces in a smoother fashion than rectangular 
or square patches. They model currents in a piece-wise linear manner and charge density 
in a pulse doublet manner. They do not introduce extraneous line charges since the current 
from one triangle must flow completely to the second. 


The fundamental vector direction of these 
basis functions is from the centroid of one 
triangle to the centroid of the adjacent 
triangle forming the couple. This leads to a 
ZIG ZAG representation of current direction 
on a body. Consider the plate geometry 
shown in figure 12-1 modeled with triangle 
couples. The triangles are shown with 
dashed lines while the current basis 
functions are shown in solid arrows from 
centroid to centroid. We see the inherent 



Figure 12-1 Basis function currents are from 
triangle centroid to centroid (zig-zag) 
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change of vector direction from one couple to the next, interior to the plate and along its 
edges. 

Also note that the basis function direction on each corner tip is different for adjacent 
corners and is a function of body symmetry. Referring to Figure 12-2, the plate with no 
symmetry has all comer couple directions parallel which makes the basis function direction 
parallel or perpendicular to the diagonal for adjacent tips. For the plate with symmetry the 
basis function vector directions are no longer all parallel and is different for each pair of 
tips. The symmetric plate has a basis function representation that is not symmetric. 

The ZIG ZAG nature of the triangle couple basis functions will change if one invokes model 
symmetry. Figure 12-2 shows a typical triangle modeling scheme with and without symmetry. 
Note that the comer tip couple vector direction is different for each of these plates. 

This ZIG ZAG nature of the vector direction is inherent in the triangle basis functions 
introduced by [1]. 



No Symmetry 


Mirror Symmetry 



Figure 12-2 Geometry symmetry alters basis vector orientation 
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Physical currents do not have this ZIG ZAG characteristic. Triangle couple representation 
of physical surface currents introduces a modeling artifact. The consequence of the ZIG-ZAG 
representation is usually minor, yet, on occasion it will introduce minor asymmetric results 
in radiation patterns. This often times is seen in edge wave scattering from a square plate 
where the computed results do not have complete symmetry as one would expect This is 
due to the tip basis functions not being the same direction on all four comers. 
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SECTION 13 


EXAMPLE RESULTS 


This section presents sample MOM3D results for backscatter RCS, bistatic RCS, current 
distributions, downrange images, down/cross range images, near zone field maps, and 
antenna gain patterns. These results were computed on a personal computer with a 486 
processor operating at 33 MHz. 

Coordinate System: The coordinate system is the standard spherical system with the 
azimuth angle <p measured from the X axis and the elevation angle measured up from the 
X-Y plane, Figure 13-1. The polar angle 0 is the 90° complement of elevation and is 
measured down from the Z axis. The two 
fundamental polarization unit vectors are 
u* and u f that are in the "direction" of the 
corresponding angle 0 and <p. The xyz 
rectangular representation of these vectors 
become functions of (0,<p). 


Model Geometry: The model geometry was 
a 5.5 inch square plate at a frequency of 3 
GHz. The side length and area were 1.4 A. 
and 1.95 A. 2 respectively. The plate was 
centered in the X-Y plane and was 
modeled with 200 triangles with mirror symmetry about the Y axis. Figure 13-2. This 
resulted in 280 unknown current couples. The symmetric geometry resulted in an "even" 
matrix rank of 145 and an odd matrix rank of 135. There were 10 current couples on the 
spine. The normalized sample density was 143 current couples per square wavelength, which 
is over-sampled. The plate was a perfect conductor. 



Figure 13-1 MOM3D coordinate system is 
the standard spherical system 
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Figure 13-2 PLT-55.GEO Mesh Model, Symmetric about Y = 0 (X-Z) plane. 


Backscatter RCS: The backscatter RCS is shown in Figure 13-3. An elevation cut is shown 
for both 6 (solid line) and <p (dashed line) polarizations. We see that the RCS is the same 
for both polarizations when perpendicular to the plate, Le., the "specular flash" where the 
angle of reflection is equal to the angle of incidence. This flash has a value of -4 dBsm that 
compares to the high frequency physical optics value of = 4* A 2 / A. 2 = -3.2 dBsm. Thus 
even for this electrically small plate we obtain optics results for the specular return. For 0 
polarization at 0° incidence the electric field is perpendicular to the plate so that no currents 
are induced (no tangential E field), and therefore no scattered electric field. 

The 0 polarization lobe at 40° is the surface traveling wave (TW) lobe. The location 
compares well to the hip pocket standard location given by 8 = 49 ( A / L ) 0J = 41°. The 
TW magnitude is -14 dBsm, which compares to the hip pocket estimate that the TW return 
be less than 3 A. 2 * -15 dBsm. 
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(a) Elevation Cut 


(b) Azimuth Cut 


Figure 13-3 Backscatter RCS 

For <p polarization at 0°, the RCS is due to edge currents, - 22 dBsm, which corresponds to 
the standard edge optical return of = L 2 / it = -22 dBsm. 

The azimuth cut, Figure 13-3(b), is for <p polarization and is typically called the edge wave 
cut since the RCS is dominated by specular edges and edge traveling waves. The edge 
specular is the same as the previous edge value, -22 dBsm. The edge traveling wave here 
is not particularly dominant The slight asymmetry in this result is due to the ZIG-ZAG 
basis function representation and the effect of model symmetry that causes asymmetric basis 
functions (see Section 12). 


Bistatic RCS: Bistatic RCS is shown in Figure 13-4 for 6 and <p polarizations (solid and 
dashed lines respectively) for an elevation cut In (a) the plate is illuminated at 90°. The 
backscatter RCS at 90° is -4 dBsm. This is the optics specular flash and corresponds to 
= -3 dBsm. The forward scattered field at 270° has the same amplitude as the backscatter 
value. This forward scattered field, when added to the incident field (in a vector phasor 
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Q 


0 



(a) 90 deg incidence (b) 45 deg incidence 


Figure 13-4 Bistatic RCS 

manor) forms the shadow behind the plate. At edge viewing angles of 0° and 180° there is 
no scattered field for 6 polarization (perpendicular to plate) but there is a finite <p polarized 
scattered field due to edge induced currents. 

The bistatic result for 45° incidence is shown in (b). The specular flash here, where angle 
of refection is equal to angle of incidence, is at 135°. The specular value is -7 dBsm and 
compares to the physical optics value of = cos 2 « 4* A 2 / A. 2 = -6 dBsm that is down by 
3 dB over the 90° specular result due to the cos 2 (45°) = 0.5 factor. The forward scatter lobe 
at 225°, (= 45° + 180°), forms the shadow behind the plate when added to the incident field. 
For edge viewing angles, as before, there is a <p component of scattered field, but no 6 
component 
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Square 5.5 Inch Plate; 200 Triangles 

CURRENT M0M30 COMPUTATION 



Figure 13-5 Current examples 


Current Computation: Plate currents are shown in Figure 13-5 for the RMS (root mean 
square) average for four excitation angles. In (a), the excitation is perpendicular to the plate 
for E* polarization. We see that the current peaks along the pair of edges parallel to the 
incident field and tends toward zero for the edges perpendicular to the incident field. 
Currents in the center of the plate correspond to physical optics levels. We note that the 
edge parallel currents show a 3 peak standing wave pattern corresponding to energy flashing 
between the vertices (recall that the plate edge length was approximately 1.4 A. that would 
create 3 standing wave peaks, one every half wavelength). The current along the edge 
perpendicular to the incident E field, while tending to zero, show an oscillatory pattern. 
Recall that triangle currents are the vector sum of the common (couple) edge currents. 
Thus referring to Figure 13-2, we see that edge triangles alternate between having two or 
three common edges. Those with three will have higher currents. 

Figure 13-5 (b) shows the RMS induced currents for the surface traveling wave excitation 
case, E* polarization at 20° elevation and 0° azimuth (looking perpendicular to edges). The 
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results show low levels of currents near the front edge that build in intensity going toward 
the rear edge. A three peak current standing wave pattern is evident, one eveiy XU. This 
is caused by the forward current phasor adding with the current reflected by the rear edge. 
We also see two pronounced edge current peaks along each plate side, again showing a 
standing wave pattern due to forward and reflected edge currents. 

Figure 13-5 (c) shows the RMS currents for parallel edge on excitation. The incident field 
is 9 polarized. We see a dominant leading edge current at the plate front edge. The side 
edges show the edge standing wave caused by the forward and reflected edge energy. The 
currents in the plate center tend to be small. The rear edge current show a small peak, 
probably due to diffraction energy flashing between the two back vertices. 

Figure 13-5 (d) shows the RMS currents for the edge wave excitation at 45° azimuth for <p 
polarization. We see two symmetric edge standing wave currents, again with a peak every 
XU. These build in amplitude as the wave travels toward the aft vertices where they diffract 
energy toward the front and rear comer vertices. The rear comer vertex current also peaks. 
Currents in the center of the plate are small since this is an edge wave dominated excitation. 

Downrange Images: Downrange images are shown in Figure 13-6. The images presented 
here correspond to co-polarized backscatter energy (recall that the method described in 
Section 5 also could compute a bistatic and/or a cross polarized image). 

Case (a) is for illumination perpendicular to the plate such that the backscatter is the 
specular flash, -4 dBsm. This orientation shows all energy originating from a single 
downrange location. 

Case (b) is for the traveling wave excitation, E* polarization at 20 ° elevation. We clearly 
see the return as originating from the aft edge of the plate. This is due to the reflection of 
the "forward” currents at the aft discontinuity. When these surface currents reflect and 
travel toward the plate front edge they reduce in intensity as they radiate energy back 
toward the illuminating source that causes the image level to decay going toward the front 
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(a) E-th, (Az,EO - 0,90 



Figure 13-6 Down range image of plate 
edge of the plate. 



(b) E-th, (Az,EI)- 0,20 
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(d) E-phi, (Az,E0- 45,0 


Case (c) is for E* front edge parallel polarization for 0° azimuth and elevation. The image 
shows the scattered energy as originating from the front edge currents and corresponds to 
the physical optics edge value of = l 2 / it = - 22 dBsm. We also see a return from the 
aft end of the plate due to reflected edge wave energy. 


Case (d) is for E* polarization at 45° azimuth that is the edge wave illumination case. We 
see a peak return corresponding to edge wave reflected energy from the comer vertices and 
a smaller return from the front vertex and the extreme aft vertex due to multiple diffraction. 


Cross Range Images: Two dimensional cross range images are shown in Figure 13-7. The 
5.5" plate dimension corresponds to ± 0.07 m for cross and downrange. These images are 
co-polarized in the backscatter direction. The cross range direction, Figure 5-1, is in the <p 
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direction that corresponds to a "rotation" about the Z axis. The images are shown as 
contour plots using PLOT 88 FORTRAN graphics. The contour lines are spaced in 2 dB 
increments. The cases presented here are the same illumination excitation as shown for the 
downrange examples, Figure 13-6, and the induced current plots, Figure 13-5. 

Case (a) is for E* illumination perpendicular to plate, (Az,El) * 0", 90°. The backscatter 
return here is the "specular” flash where the plate currents have mostly the same phase and 
when radiating, phasor add to form the "flash." The two dimensional image shows radiation 
as originating at a constant downrange distance distributed over the cross range dimension 
of die plate. Recall that the specular flash had a value of -4 dBsm while the image shows 
a maximum contour of between -14 and -12 dBsm. Two dimensional images thus show 
distributed levels. 

Case (b) is for the surface traveling wave illumination, E* polarization at 20° elevation and 
0° azimuth (perpendicular to edges). As expected, the image shows the major scattering as 
originating from the surface traveling wave reflecting from the aft plate edge. Note that the 
intensity then falls off as this reflected energy decays as it flows back toward the front edge 
(Figure 13-6b). The intensity is slightly higher along the downrange edges. 

Case (c) is for parallel edge illumination, E* polarization at 0° azimuth and elevation. The 
image clearly shows the leading edge as the dominant scattering source due to leading edge 
induced currents. The two rear vertices are lesser scattering centers due to energy flowing 
down each downrange edge. The source at the center of the rear edge is probably due to 
multiple diffraction from the rear vertices. 

Case (d) is the edge wave illumination example, E* polarization at 45° azimuth and 0° 
elevation. The major scattering centers are the two vertices that reflect the edge waves. 
The far rear vertex also is a major scattering center due to multiple vertex diffraction of the 
edge traveling waves. The forward vertex is not a significant scattering center. Compare to 
the downrange result, Figure 13-6d. 
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Figure 13-7 Two dimensional range / cross range images 
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Near Field Examples: Near field 

computations for E* 0 * - and E“* are shown 
in Figures 13-9 and 13-10. The plate 
excitation was E* perpendicular to the 
plate, 90° elevation, 0° azimuth. The near 
zone fields were computed over a 6 A. 
square grid centered in the X-Z plane, 
Figure 13-8, with computation points 


computed, E total and E scattered. For 

each, the root mean square (RMS) scalar time average was computed along with the two 
principal time values, one for cot = 0° (the REAL part), and the second for cot = -90° (the 
IMAG part). The results displayed are for contours of intensity, not vector direction. 



E Total Field: The total electric field. Figure 13-9, is the vector sum of the incident and 
scattered field, E To * = E“ c + E“*. The incident field for this case is a plane wave 
propagating down the Z axis. 

Case (a) is the RMS average. Beneath the plate the incident and scattered fields are out 
of phase to form a shadow region. Above the plate is a standing wave pattern every A./2. 
The peaks of the standing wave occur when the scattered and incident fields are in phase 
while the nulls occur when the two fields are out of phase. The first peak occurs at a height 
above the plate of A./4. 


Case (b) is the IMAG value of the field. This is an instantaneous snap shot in time 
corresponding to cot = - 90°. The alternating intensity of the incident plane wave is clearly 
seen. Above the plate the reflected (scattered) field causes an interference pattern (standing 
wave) with the incident field. Beneath the plate the interference pattern forms the shadow 
region. 
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Figure 13-9 Total field 

Case (c) is the REAL value of the field. This corresponds to a snap-shot in time of cot = 
0°. The incident wave for this case has moved A/4, i.e., 90°, relative to case (b). Thus the 
peaks of the incident field for (c) are located where the nulls in (b) are located. 

E Scattered Field: The scattered field E”* is shown in Figure 13-10. The equivalent far 
field RCS pattern for this case is the bistatic result shown in Figure 13-4a. The backscatter 
direction (reflected) is above the plate while the forward scattered field, which forms the 
shadow region beneath the plate by phasor subtraction with the incident field, is beneath the 
plate. The back and forward scattered fields are symmetric above and below the plate. 

Case (a) shows the RMS average intensity contours. The four lobes at 45° intervals are also 
seen in the bistatic result, Figure 13-4a, at 45, 135, 225, and 315 degrees. The field is most 
intense near the plate. 
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E Scattered Reid 

(a) RMS Average 

(b) Imag ( at - -90 s ) 

(c) Real (ot-O* ) 


Case (b) and (c) show the two instantaneous snap-shots of the scattered field. The two are 
separated by 90° phase, i.e., JJ4, such that in (b) the peaks occur where the nulls occur in 
(c). Far from the plate the field lines close back on themselves, i.e n the field is solenoidal. 
Near the plate the field lines terminate on the induced charges. 


Antenna Example: MOM3D also may 
be used for antenna analysis. This 
example shows the computed gain 
pattern for a half wave dipole modeled 
as a thin strip, Figure 13-11. The dipole 
is centered along the Y axis in the X-Y 
plane. The center port (couple) is the 
feed location, 8. The computed gain 
pattern. Figure 13-12, shows the typical 
half wave dipole figure eight pattern 
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Figure 13-11 Half wave dipole strip antenna 
geometry 
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Figure 13-12 Half wave dipole antenna gain pattern 

with peak gain of 2.1 dbi. The computed input impedance was Z = 72 + j 6 ohms that 
corresponds to a resonant half wave dipole operated slightly above resonance since the 
reactive component is inductive. 
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APPENDIX A 


NUMERICAL and ANALYTICAL MATRIX ELEMENT INTEGRATION 

The approach and theory used to compute the Green’s function integrals in a less 
approximate or more exact manner than the centroid approach is presented. MOM3D uses 
die centroid approximation for all triangle interactions greater than approximately ten 
triangle spacings apart For closer interaction distances, as specified by the user, one may 
choose from three approaches for computing matrix elements: the centroid approximation; 
or a combination analytical / numerical approach where the numerical integration may use 
a single point (centroid) or may use three points. 

Theory is presented for: numerical quadrature formulas for surfaces; singularly removal; and 
an analytical formula for the 1/R integration (point in space integrated over a triangular 
area). 


Overview: 


The initial 



development of MOM3D used a 
centroid approximation when 
computing triangle to triangle 
matrix element interactions, 

Figure A-l. This results in a 
considerable reduction in i n i . i rc . r ( - m i 

complexity and in the computer 

* * r Distance between source and field points written in terms of the 

. • . , . - . centroid distance. 

time required to fill the im- 
pedance matrix. This technique 

is reasonable when the mesh _ 

Figure A-l Interaction distance R 

triangles are greater than one 

characteristic mesh length apart Early work showed that the approach was still viable even 
for touching co-planar triangles. But, the centroid approximation becomes less reasonable 
when two triangles become very close to each other, such as when one triangle bends back 


Distance between source and field points written in terms of the 
centroid distance. 


Figure A-l Interaction distance R 
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on the second, as in modeling a small angle wedge, or when modeling two parallel surfaces 
with separation less than a characteristic dimension L. The fundamental difficulty is with 
the 1/R integrand. When R is greater than L, the centroid approximation is reasonable, but 
when the distance between triangles is less than L, then the centroid approximation over 
estimates the 1/R interaction. Early studies showed that the centroid approximation was a 
good approximation when triangle spacing was greater than one characteristic dimension L 
apart, i.e., depends on mesh size and not on wavelength. Further, the centroid 
approximation was adequate for two triangles with a common edge whenever the included 
angle between was greater than 50 degrees, i.e., valid for planar surfaces and those with not 
too rapidly varying curvature. 


The basic current unknown is two adjacent 
triangles with a common side, which forms 
a couple. The matrix element interaction 
is between pairs of couples, Figure A-2, and 
involves interaction between four triangles. 
The impedance matrix element between the 
ith and jth couples is a sum of four terms 
accounting for the interaction between the 
four triangles making up the two couples: 




Figure A-2 Matrix elements involve 4 
triangle interactions 


z, • <w„n})-. 


EE hb s , s » 


; 0) G w> 
7 w« 


(A-l) 


A-2 



This expression is exact in the sense that no approximations have been made other than the 
fundamental choice of basis and weight functions. The terms G (1) and G (0) involve 
integration of the Green's function between pairs of triangles, 
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These are four fold integrals that express the electrical interaction of one surface patch with 
another. 

The centroid approximation for G° and G 1 is: 


e -jtr 


4* \r\ 


where f c * *i~T) » and 
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where ^ is the distance between the centroid of the ith and jth triangles. This approach is 
computationally simple and fast, which is important when solving large geometries. This 
approximation begins to fail when the triangle interaction distance becomes less than a 
characteristic dimension such as for very close parallel surfaces or for knife edges where one 
triangle folds back sharply on another. 
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Numerical Triangle Integration Formulas: An integration scheme similar to one 

dimensional Gaussian quadrature, but for triangular surfaces and tetrahedron volumes is 
discussed by Hammer, Marlowe, and Stroud [Al]. They derive general surface and volume 
numerical integration formulas of the form 


//«) dV - £ aj/Uj) 
7-1 
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where the numbers are constants (weights) (j=l to k), ^ are points in the domain of f, 
and R is the bounded closure of an open set in E„ where n is the dimensionality of the 
space. Their goal was to obtain numerical integration formulas for the n-simplex to hold 
exactly for polynomial functions f of at most degree k. 

Hammer et al. [Al] give weight sets, aj, and function evaluation points, £j, for 1, 3, 4, and 
7 function evaluation points over a triangle. The three point formula is said to hold exactly 
for quadratic functions, the four point formula for cubic functions, and the seven point 
formula for quintic functions. 

It is interesting that reference [Al] specifies the centroid as the first affine invariant formula 
for the triangle as the sole evaluation point with the weight equal to the area. This of 
course is our centroid approach. 

Application of reference [Al] results requires us to represent the function evaluation points 
as vectors. Writing the vertices of a triangle as V„ V* V 3 and the centroid as C = 1/3 £ Vj, 
the function evaluation points r r (inside the triangle) are specified as 


r? t + (1 - r)€ where i - 1, 2, 3. 
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where r is a constant specified by [Al] for the various orders of numerical integration. The 
function evaluation points inside the triangle can be rewritten in terms of the triangle 
centroid C as 


f f - C - r(C - fy 


(A-6) 


where we recognize that the vector ( C - V s ) is from the i th triangle vertex to the centroid, 
Figure A-3. 



C: Triangle Centroid 

Vi : i th Triangle Vertex 

C - VI : Vector from vertex to 
centroid 


VZ 



vz 



VI 


r = C - r ( C • Vi ) 


Figure A-3 Triangle function evaluation points 

Three point surface integration is exact for quadratic functions (e.g., 2 terms of a Taylor 
series expansion) for r = ± 1/2 with weights W - 1/3 A where A is the triangle area. For 
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r = + 1/2, the evaluation points are the distinct trisection of the median chords and for r 
= - 1/2, the evaluation points are the midpoints of the sides. 

Four point surface integration is exact for a cubic polynomial for r = 2/5 with weights 25/48 
A for each of the three median chords and the centroid, r * 0, with weight (-9/16) A. 

Seven point surface integration is exact for a quintic polynomial for function evaluations at 
the centroid with weight (9/40) A and by two points on each median chord given by 

r m ( 1 + v!5 ) / 7 with weight ( 155 - V15 ) / 1200 A 
r * ( 1 - V15 ) / 7 with weight ( 155 + V15 ) / 1200 A 

or approximately 

r * +0.6961405 with weight 0.1259392 A 

r = -0.4104262 with weight 0.1323942 A 

The astute reader will recognize that the sum of the weights for the 3, 4, or 7 point formulas 
is just the area of the triangle A. 

Figure A-3 shows the evaluation points within a triangle for the 3, 4, and 7 point integration 
formulas. 

Application to the double surface integral of MOM3D requires triangle to triangle 
integration, Le., a double sum. Thus the 3 point formulas will require 9 function evaluations, 
the 4 point formula 16 evaluations and the 7 point formula 49 evaluations. 

Analytical Point to Plane Line Integral Approach: Another method of computing surface 
integrals is to use Stoke's theorem to reduce an area integration to a peripheral line 
integral. If we let f represent the function to be integrated over a surface then: 
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where the function f must be represented by the curl of the vector n (vector potential). 
Stoke's theorem transforms the surface area integration to a closed perimeter line integral. 

This technique has been discussed by Baker and Copson [A2] in terms of diffraction of light 
by apertures in terms of the Maggi transformation, and by Gordon [A3] who applied the 
technique to Physical Optics scattering by replacing a surface integration with the equivalent 
peripheral line integral. 

The resulting line integral can then be integrated analytically, as done by Gordon [A3] for 
PO scattering by triangles or numerically, e.g., by Gaussian quadrature. 

The key to applying this technique is to find the vector function n for a given integrand f. 
This is not trivial. We would need a vector n corresponding to f = e - * r / r. An alternative 
approach is to remove the 1/r singularity and develop a vector n corresponding to f = 1/r. 

Singularity Removal and 1 R Formulas: The main difficulty in the centroid approximation 
is the dominance of the 1/R singularity of the Green's function. Several authors have 
subtracted out this singularity and split the Green's function integration into two terms, one 
with only the 1/R singular part and the other with the singularity removed. The benefit of 
this approach is that the term with the singularity removed is a very well behaved function 
with a defined limiting value when R goes to zero. This well behaved function can be 
numerically integrated while the singular part is dealt with separately. The following 
approach follows that by Rao [A4] who developed an analytical formula for triangles for the 
1/R singular integration. 

The starting point is to subtract and add 1/R to the Green's function resulting in 
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We will loosely call the first non singular term the phase. We note that this term has the 
finite value of -jk when R = 0 as can readily be determined from its series expansion. Thus 
the phase term is well behaved at R = 0 and can readily be integrated using any desired 
numerical scheme. 

The goal now is to integrate the singular 1/R term. Following [A4], a unique cylindrical 


Field point r f 



V 


HI 

Figure A-4 Coordinate system for 1/R triangle integration 

coordinate system (p,<p,w) is introduced. Figure A-4, where the polar axis w is through the 
field point r f and is parallel to n, the unit normal to the plane of the triangle. The w = 0 
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plane is coincident with the plane of the triangle. The field point r f isatw=±d(d 
assumed positive) along the w axis. The singular integral is written, using wn=l, as 
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where p is the radial distance to r' in the triangle. This area integral is then converted to 
line integral using Stake's theorem, 


f Vxfl 'A dS * fi-dl (A-10) 


where the line integral is completely around the triangle. The vector n is chosen to have 
only a <p component (in the plane of the triangle), thus 
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Integrating, [A4] obtains 
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where 9 is the unit vector in the 9 direction (right hand rule applies). This function is 
singular when p = 0 is included in the integration, which happens whenever the field point 
r f is on or inside the parallel projection of the triangle parallel to n or w. 


Singularity when FMd Point is inside projection of triangle 


A Singularity wh«n Field Point is on edge projection of triangis 

■ 

A — 



Figure A-5 Singularity removal depends on projection of field point r f 
The singularity is excluded by not including it in the integration, Figure A-5, 



(A-13) 


where 
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The singularity integration, where p = € => 0, is 


J, - lim f + d 2 d$ - yd 


(A-15) 


where y is the amount of polar angle required to integrate around the singularity. 

If the field point r r lies outside the projection of the triangle, no singularity exists and y = 
0. If r r lies interior to the projection of the triangle, y = 2n; if r f lies on one line of the 
triangle, y = ic; and if r f lies on two lines of the projected triangle (i.e., a vertex), then y = 
the included angle of the two lines forming the vertex, Figure A-5. 

We are now left with performing the integration of 
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where we note that $-dl = p d<J>. This line integral could be integrated numerically, e.g., 
Gaussian quadrature. However, Rao [A4] has analytically performed the integration for 
the case of a triangle. The integration is along each side of the triangle from <pj to <p i+1> 
which are the polar angles corresponding to the end points (triangle vertices) of each line 
segment The interested reader is directed to [A4] for the details. 
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Before presenting the result, we define the relevant parameters according to Figure A-4, 
which differs slightly from that used in [A4]. Let be the out of plane vector from each 
vertex to the field point r f . Let ft be the in (triangle) plane distance from each vertex to the 
w axis (p=0). Let p 0 be the in plane perpendicular distance of closest approach from the 
w axis to the line segment i to i+1. Let r a be the location on each line segment where p„ 
intercepts the line segment i to i+1 (Note that this point can be outside the triangle on the 
line extension). Let a ; and ft be the positive scalar lengths of the corresponding vector. Let 
Ij, U; be the unit vector along and in plane perpendicular to the ith edge. Referring to 
Figure A-4, the following definitions should be evident: 
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The final analytical formula for the singular function 1/R over a triangle is: 
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where the singularity removal is 


y ■ 0, 2ic, ic, or a 


(A-19) 


This analytical result enables us to compute the singular portion of the integral. 

When the field point r ( is far from the source triangle, this formula should not be used since 
in contains the difference of large quantities. However, at large distances, we will simply use 
the centroid approach. 

Combined Numerical and Analytical Result: MOM3D has used these results by using the 
1 and 3 point numerical integration formulas with the singular 1/R analytical formula to 
compute G°g. Higher order numerical integration could have been implemented at the 
expense of additional matrix fill time. It was decided to start simple and implement better 
formulas only if required. 
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The single point formula uses the triangle centroid as the numerical integration function 
evaluation point, as specified by [Al], for the non singular "phase" term and the 1/R 
analytical formula: 



- 1 + *,(?) ♦ T,(?) 
Wl 2 
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where . 

The three point formula uses the 3 interior points specified by [Al] for integration. The 
"phase" term is a double numerical integration while the singular term is a single numerical 
integration that is averaged to enforce symmetry: 
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